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Abstract 

Methods  to  estimate  the  optical  transfer  function  (OTF)  of 
the  degradation  and  the  noise-to-signal  ratio  in  a  degraded  im^ge 
are  presented  which  use  no  prior  information  about  the  degradation 
and  the  noise.  Two  types  of  degradations  are  considered: 
uniform  linear  camera  motion  and  out-of-focus  camera  with  circular 
aperture.  The  OTF  of  the  out-of-focus  camera  is  estimated  by 
inspection  of  the  spectra  of  the  degraded  image  and  from  average 
over  the  spectra  of  subsections  of  the  degraded  image.  A  line  es¬ 
timation  method  is  described  to  estimate  the  degrading  OTF  of  the 
linear  camera  motion  in  one  direction.  An  additional  method  called 
the  "logarithmic  estimate"  is  also  developed.  Frequency  and  space 
domain  estimation  methods  are  developed  for  the  noise-to-signal  ratio 
for  the  case  of  additive  gaussian  noise.  Inverse,  Wiener,  and  power 
spectrum  equalization  filters  based  on  these  estimates  are 
implemented  on  a  digital  computer  to  restore  a  variety  of 
degraded  images  and  examples  of  blind  restoration  of  these  degraded 
images  are  presented. 


PARAMETER  ESTIMATION  FOR  THE 
BLIND  RESTORATION  OF 
BLURRED  IMAGERY 


I.  Introduction 


Background 

One  of  the  major  applications  of  digital  image  processing  is 
image  restoration.  Image  restoration  is  an  estimation  process  that 
attempts  to  recover  an  original  (ideal)  image  from  a  degraded  image. 
These  degradations  are  incurred  while  the  image  was  being  acquired 
and  may  include  the  blurring  introduced  by  optical  systems,  atmospheric 
turbulence,  camera/object  motion,  as  well  as  distortion  and  noise  due 
to  electronic  and  photometric  sources. 

The  overall  performance  of  image  restoration  process  greatly 
benefit-''  from  modelling  of  the  transfer  function  of  the  degradation 
and  characterization  of  the  noise  in  the  degraded  image.  If  a  prior 
knowledge  of  the  image  and  the  imaging  system  is  available,  it  can  be 
used  to  increase  the  performance  of  the  restoration  system.  If  little 
or  nothing  is  known  about  the  image,  one  can  attempt  to  model  and 
characterize  the  sources  of  degradation  (blurring  and  noise)  and 
subsequently  remove  or  reduce  their  effects.  When  a  prior  knowledge 
about  the  image  isnot  available,  information  about  the  degradation  must 
be  extracted  from  the  observed  (degraded)  image.  This  task  is  called 
blind  image  restoration. 


The  Foreign  Technology  Division  (FTD)  at  Wright-Patterson  AFB 
has  needed  a  blind  image  restoration  method  to  deblur  photographs  taken 
under  controlled  conditions.  Most  of  the  restoration  filters  such  as 
Wiener  filter,  inverse  filter,  power  spectrum  equalization  filter 
require  the  degrading  optical  transfer  function  and/or  the  noise-to- 
signal  ratio  as  their  parameters.  In  order  to  restore  a  degraded  image, 
one  must  know  the  degrading  OTF  (optical  transfer  function)  and  the  NSR 
(noise-to-signal  power  ratio).  Since  no  prior  information  about  these 
two  parameters  is  available,  one  has  to  estimate  them  from  the  degraded 
image  itself.  Once  the  degrading  OTF  and  the  NSR  are  estimated,  any  of 
the  restoration  filters  (methods)  can  be  employed  to  restore  the  degraded 
image.  Therefore,  the  problem  of  the  blind  restoration  reduces  to 
estimation  of  the  degrading  OTF  and  the  NSR. 

Statement  of  the  Problem 

The  problem  that  is  addressed  in  this  thesis  may  be  stated  as 
follows : 

Given  a  degraded  image  with  no  information 
available  concerning  the  degradation  function 
and  the  noise,  develop  methods  to  estimate  the 
degrading  OTF  and  the  NSR  and  use  them  to 
restore  the  degraded  image. 


Scope 

This  thesis  is  primarily  concerned  with  the  estimation  of  the 
degrading  OTF  and  the  NSR.  Inverse  filter,  Wiener  filter,  and  power 


2 


spectrum  equalization  filter  are  presented  as  restoration  methods. 

Other  restoration  techniques  and  the  investigation  of  their  performance 
with  the  estimated  degrading  OTF  and  the  NSR  are  considered  out  of  the 
scope  of  this  thesis.  Two  specific  cases  of  the  degradation  will  be 
considered:  uniform  camera  motion  and  out-of-focus  lens  system  with 
circular  aperture.  Restoration  for  other  forms  of  degradations,  such 
as  blurring  due  to  atmospheric  turbulence,  camera  vibration,  geometric 
distortion,  and  system  nonlinearities  will  not  be  addressed.  The 
primary  objectives  of  this  effort  can  be  stated  as: 

1.  Develop  a  method  to  estimate  the  degrading  OTF. 

2.  Develop  a  method  to  estimate  the  noise-to-signal 
power  ratio. 

3.  Implement  the  results  by  a  general  purpose 
restoration  filter  on  a  digital  computer. 


Assumptions 

The  imaging  system  and  the  degradation  system  will  be  assumed 
as  linear  space-invariant  systems.  We  will  assume  that  the 
original  image  and  the  noise  are  spatial  stochastic  random  processes 
and  that  they  are  spatially  stationary  and  statistically  independent. 
The  noise  process  is  assumed  to  be  additive  gaussian  process  with  mean 


Approach 

Presentation  of  this  research  will  proceed  in  Chapter  2  with 
a  discussion  of  general  theory  concerning  the  linear  space  invariant 
system,  image  formation  system,  the  degrading  optical  transfer  functions, 
noise  process  and  image  restoration  methods.  This  is  followed  in 
Chapter  3  with  discussions  of  the  methods  to  estimate  the  degrading  OTF 
for  the  cases  of  uniform  linear  camera  motion  blur  and  out-of-focus 
camera  blur  and  the  noise-to-signal  power  ratio.  In  Chapter  4,  we 
discuss  the  implementation  of  these  estimation  methods  and  their 
performance  in  actual  degraded  images.  In  Chapter  5,  we  summarize 
our  findings  and  recommend  topics  for  continuing  research. 


Materials  and  Equipment 

All  thesis  research  was  accomplished  in  the  Signal  Processing 
Laboratory  at  AFIT  utilizing  the  Data  General  Eclipse  S/250  and 
Nova  computers  and  the  Octek  Image  Analyzer.  The  digital  images 
used  in  this  research  were  256x256  monochrome  images  with  four  bits  per 
pixel.  The  images  were  acquired  by  a  variety  of  techniques,  including 
optical  scanning/digitization  of  photographic  film,  frame  grabbing  from 
a  vidicon,  and  computer  generation.  Programming  was  accomplished 
using  FORTRAN  IV  and  FORTRAN  5.  Existing  software  was  used  for  image 
acquisition  and  display,  two  dimensional  Fourier  transformation  and 
inversion,  and  image  file  input/output. 


II.  Image  Restoration  Theor 


Linear  Space  Invariant  Systems 

We  can  think  of  a  system  as  a  mapping  of  a  set  of  input 
functions  into  a  set  of  output  functions.  For  the  case  of 
electrical  networks,  the  input  and  output  functions  are  real 
functions  (voltages  or  currents)  of  a  one  dimensional  indepen¬ 
dent  variable  (time):  for  the  case  of  imaging  systems,  the 
inputs  and  outputs  can  be  real-valued  functions  (intensity)  or 
complex-valued  functions  of  a  two-dimensional  independent 
variable  (space). 

A  single-input  single-output  two  dimensional  system  can  be 
represented  by  a  mathematical  operator  "P",  which  operates  on  an 
input  function  f(x,y)  to  produce  an  output  function; 
g(x,y),  i.e., 

g(x, y)  =  P  { f (x, y)}  (1) 

The  system  is  linear,  if  the  superposition  property  is 
obeyed  for  all  input  functions  f^(x,y),  and  f2(x,y),  and  all  com¬ 
plex  constants  "a"  and  "b"; 

P{  af  j  (x ,  y )  +  bf2(x,y)}  =  aPtf^x.y)}  +  bP{f2(x,y)}  (2) 

The  great  importance  of  linearity  is  the  ability  to  express  the 
response  to  a  complicated  input  by  first  decomposing  this  input 
into  a  linear  combination  of  elementary  functions  and  then  taking 
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the  same  combination  of  these  elementary  responses.  Such  a 
decomposition  may  be  obtained  by  the  sifting  property  of  the  delta 
function  which  states  that 

CO  CO 

f (x , y )  =  S  I  f(Y,8)  <5(x-Y,  y-6)  dY  dB  (3) 

-00  -CO 

f(Y.  3)  is  considered  as  weighting  factor  applied  to  6(x-y,  y-8). 
Combining  Equations  (1)  and  (3),  we  obtain 

00  00 

g(x, y)  =  P{  I  f  f(Y,8)  6(x-y,  y-8)  dy  d6)  (4) 

_°°  -00 

If  the  system  is  a  linear  system,  then  applying  linearity  property 
to  Equation  (4),  we  can  obtain 


g(x,y)  =  f  f  f(Y,6)  P{ 6(x'Y»  y-8)}  dY  dB  (5) 

—CO  _oo 

Now,  let  the  h(x,y;  Y,8)  be  the  response  of  the  system  at  point  (x,y) 
of  the  output  space  to  a  delta  function  input  at  coordinates  (y>8)  of  the 
input  space,  thus, 

h(x,y;  y,8)  =  P{<5(x~y,  y-8)}  (6) 

This  function  is  called  the  point  spread  function  (PSF)  of  the 
system.  The  system  input  and  output  now  can  be  related  to  each 
other  by  the  following  equation 


OO  00 


g(x,y)  =  /  /  f(Y,8)  h(x,y;  Y,8)  dY  d8 


(7) 
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This  expression  is  known  as  the  superposition  integral. 

A  system  is  said  to  be  space-invariant  (shift-invariant  or 
isoplanadc ) ,  if  the  behavior  of  the  system  is  not  a  function  of 
the  independent  variable  (space).  Thus  if  the  linear  system  is  also 
space-invariant,  then  the  impulse  response  of  the  system  will  be 

h(x,y;  y,8)  =  h(x-y,  y-S)  (8) 

and,  we  can  express  the  input  output  relationship  for  a  linear  space 
invariant  system  as  follows 

OO  GO 

g(x,y)  =  f  f  f(Y,B)  h(x-Y,  y-6)  dYd8  (9) 

—  OO  — OO 

Equation  (9)  states  that  the  output  of  a  two-dimensional  linear, 
shift-invariant  system  is  the  two-dimensional  convolution  of  the 
input  function  with  the  point  spread  function  of  the  system. 

Imaging  System 

An  imaging  system  consisting  of  lenses,  prisms,  mirrors  and  so  on 
can  be  considered  as  a  black  box  which  provides  a  transformation  or 
mapping  of  input  light  distribution  or  input  energy  radiated  from 
the  object  to  some  output  spatial  light  distribution  or  radiant 


energy. 


Entrance 

Pupil 


imaging  system 


Figure  1  Generalized  Imaging  Systems  (Ref.  1) 

Such  an  imaging  system  can  be  characterized  by  an  entrance  pupil, 
representing  a  finite  aperture  through  which  radiant  energy  (light) 
passes  through  to  reach  imaging  elements,  and  an  exit  pupil 
representing  again  a  finite  aperture  which  radiant  energy  must  pass 
as  it  leaves  the  imaging  elements  to  reach  the  image  plane. 

In  Figure  1,  a  point  source  in  the  object  plane  at  coordinates 
(x,  y)  of  intensity  f(x,y)  radiates  energy  toward  the  imaging  system 
entrance  pupil.  The  energy  coming  from  the  system  exist  pupil 
produces  an  intensity  distribution  g(x^,y^)  in  the  image  plane.  It  is 
assumed  that  passage  of  the  energy  between  the  entrance  and  exit  pupil 
can  be  described  by  geometrical  optics. 

In  most  image  formation  systems  the  energy  radiation  emitted 


by  an  object  arises  from  transmitted  or  reflected  light  from  an 
incoherent  light  source.  The  energy  radiation  can  often  be  regarded 


as  quasimonochromatic  in  the  sense  that  the  spectral  bandwith  of 
the  energy  radiation  detected  at  the  image  plane  is  small  with  respect 
to  the  center  of  the  radiation.  Under  these  assumptions,  the  imaging 
system  of  Figure  1  will  respond  as  a  linear  system  in  terms  of  the 
intensity  of  its  input  and  output  fields  (Ref.  1).  The  relationship 
between  the  image  intensity  and  the  object  intensity  for  the  imaging 
system  can  then  be  expressed  by  the  superposition  integral;  i.e., 

OO  00 

g(xi,yi)  =  /  /  h(xi,yi;  x.y)  f(x,y)  dx  dy  (1 

—00  — OO 

where  hCx^y^;  x,y)  is  the  image  intensity  response  of  the  system  at 
point  (x^.y^)  of  the  image  plane  to  a  point  source  of  light  at 
coordinates  (x.y)  of  the  object  plane.  If  the  intensity  impulse 
response  is  space  invariant,  the  input  output  relationship  can  be 
expressed  by  convolution  equation 


g(xi,yi)  =  /  /  h(xi-x,  yi-y)  f(x,y)  dx  dy  (1 

—  00  —  00 

If  we  take  Fourier  transform  of  the  both  sides  of  Eq.  (11) 

G(u, v)  =  H(u, v)  F(u, v)  (1 

where  H(u,v)  is  called  the  optical  transfer  function  (OTF),  and  the 
magnitude  |  H(u,v)  |  of  the  OTF  is  known  as  the  modulation  transfer 
function  of  the  system  (Ref.  1). 


Figure  2  Image  Restoration  System 


An  image  restoration  system  can  be  modeled  as  in  Figure  2. 

The  original  (ideal)  image  f(x,y)  is  degraded  by  the  operation 
h^Cx.y)  and  a  noise  n(x,y)  added  to  form  the  degraded  image  g(x,y). 

This  is  convolved  with  the  restoration  filter  point  spread  function 

A 

h^(x , y )  to  produce  the  restored  image  f(x,y). 

The  degradation  system  models  the  source  of  degradation  in  the 
image.  There  are  many  sources  of  degradation.  Some  types  do  not  involve 
blur,  while  effecting  only  the  gray  levels  of  the  individual  picture 
points.  Other  types  which  do  involve  blur  are  called  spatial  degra¬ 
dations.  Such  degradations  can  be  introduced  by  atmospheric  turbulence, 
aberrations  of  the  camera  and  the  object.  We  will  only  be  concerned 
with  the  spatial  degradations  which  are  due  to  uniform  linear  camera 
motion  and  out-of-focus  lens  system.  The  degradation  system  and 
restoration  system  can  be  considered  as  the  imaging  system  of  Figure  1. 
Therefore,  these  systems  will  respond  as  a  linear  system  in  terms  of  the 


intensity  of  their  input  and  output  fields.  If  the  degradation 
system  models  an  out-of-focus  camera,  then  this  system  will  introduce 
a  defocusing  at  the  output.  A  shift  in  the  position  of  its  input  will 
not  effect  the  defocusing  introduced  by  the  system.  The  only  effect 
caused  by  a  shift  in  the  position  of  its  input  will  be  an  equal  shift 
in  the  position  of  its  output,  i.e.,  it  will  introduce  the  same 
amount  of  defocusing.  This  may  be  an  idealization  when  applied  to 
physical  systems.  But,  we  will  assume  that  the  degradation  and 
restoration  systems  are  linear  spacp-invariant  systems.  Therefore, 
degraded  and  restored  images  can  be  expressed  as 

oo  oo 

g(x ,  y )  =  f  S  h,(x-Y,  y-6)  f(Y,S)  dYdS  +  n(x,y)  (13) 

—  00—00 

and 


*  oo  oo 

f(x,y)  =  /  /  h  (x-y,  y-6)  g(  Y,  0)  dYdS  (14) 

—  00—00 

By  taking  the  Fourier  transforms  of  both  sides  of  these  equations, 
we  obtain 

G(u, v)  =  Hd(u,v)  F(u, v)  +  N(u, v)  (15) 

A 

F(u,v)  =  Hr(u,v)  G(u , v )  (16) 

A 

where  G(u,v),  F(u,v),  F(u,v),  H^u.v),  H  (u,v),  N(u,v)  are  the  Fourier 

A 

transforms  of  g(x,y),  f(x,y),  f(x,y),  hd(x,y),  h  (x,y),  n(x,y) 
respectively.  The  function  Hd(u,v)  is  the  transfer  function  of  the 


1^.1 


system  that  transforms  ideal  (original)  image,  f(x,y),  into  the 
degraded  image  g(x,y).  The  function  H^Cu.v)  is  the  transfer  function 
of  the  restoration  system  which  makes  as  good  an  estimate  as 
possible  of  the  original  image  f(x,y). 

In  order  to  determine  such  a  transfer  function,  Hr(u,v)  we 
need  some  form  of  information  about  the  degrading  transfer  function. 
In  some  cases,  the  form  of  the  degradation  can  be  used  for 
determining  the  H^(u,v).  Therefore,  we  will,  first,  restrict  our 
attention  to  some  specific  forms  of  degradations  and  their  transfer 
functions.  We  will  later  consider  the  case  of  blind  restoration 
where  we  assume  no  prior  knowledge  of  the  degradation  and  must 
estimate  the  degrading  transfer  function  from  the  degraded  image. 


OTF  of  Some  Soecific  Deeradations 


Some  common  forms  of  degradations  can  be  characterized  by 
specific  transfer  functions.  If  we  know  the  degradation  form,  we 
can  obtain  the  parameters  of  that  degrading  transfer  function  from 
the  degraded  image,  and  use  these  parameters  to  obtain  the  same  form 
of  degrading  transfer  function.  We  will  examine  the  degrading  transfer 
functions  which  are  produced  by  linear  camera  motion,  and  out-of-focus 
camera. 


OTF  of  Linear  Camera  Motion 


A  relative  motion  between  camera  and  object  produces  a  blur 
in  the  image  in  the  direction  of  the  motion.  We  will  assume  that  the 
image  is  invariant  in  time  except  for  the  motion.  If  a  point  FL  is 
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selected  in  the  image  plane,  the  instantaneous  image  intensity 
at  this  point  will  be  denoted  by  g(K^,t).  The  total  exposure  at  any 
point  of  the  recorded  image  can  be  described  by  the  integrated 
image  intensity  at  this  position  (Ref.  2). 

T 

g(K.)  =  /  g(K.,t)  dt  (1 

1  0 

where  T  is  the  total  exposure  time.  The  medium  that  actually 
performs  the  integration  could  be  photographic  film  or  the 
phosphor  in  the  camera  sensor.  Assuming  that  energy  is  conserved 
by  the  imaging  system  at  any  instant  during  the  exposure,  the  energy 
radiated  and  collected  by  a  small  area  in  the  object  plane  around 
the  object  point  Kq  is  described  by 

g(x.,y.)  dx.  dy.  =  f(x  ,y  )  dx  dy  (1 

s  1  3 i'  i  3 \  v  o  Jo  o  Jo 

at  a  fixed  time  t.  And,  we  can  say  that 

g(K.,t)  =  f(Ko,t)  (1 

and  time  integral  can  be  written  as 
T 

g(K.)  =  /  f (K  , t)  dt  (2i 

0  ° 

This  relation  shows  that  a  recorded  image  point  g(K^)  over  the 
exposure  time  T  will  be  obtained  by  performing  a  summation  over  all 
object  points  that  are  imaged  onto  the  image  point.  If  there  were 
no  motion  of  either  object  or  camera,  this  relationship  would  become 


(2 


g(K.)  =  f(Ko,t)  T 

Since  we  have  assumed  that  object  is  a  two-dimensional  entity  and 

the  image  is  time  invariant  except  for  the  motion,  we  can  rewrite 

the  recorded  image  intensity  as 

T 

g(x-,yi)  =  /  f[x  (t),  y  ( t )  ]  dt  (2 

0  0 

Given  a  description  of  the  object  motion  R(t),  we  can  convert 
this  time  integral  into  a  positional  integral  over  an  equivalent 
stationary  object.  Let  the  path  of  the  object  motion  be  described 
in  two-dimensional  object  plane  as  a  parametric  function  of  the 
time 

R(t)  =  [xQ(t),  yQ(t)  ]  (2: 

The  elemental  length  ds  of  the  path  R(t)  is  given  by 


The  time  integral  in  equation  (22)  can  then  be  changed  into  the  lii 
integral  of  the  form 


R(T) 

g(*H -y J  =  f 


f[xo(t),  yQ(t)  ]  ds 


(2 


which  is  the  general  expression  for  motion  degradation  with 
linear  space  invariant  imaging  system.  The  space  invariant 
point  spread  function  can  be  identified  in  Equation  (25)  as 


h(x . , y . ;  x  , y  )  = 

i  3 1  o  o 


'dx^t) 
dt 


dyjt) 
dt 


(21 


Where  x  (t),  yQ(t)  are  valid  over  the  path  of  the  object 
motion,  that  is, 

xQ(0)  1  x0(t)  1  xo(T) 

(2 

y0(°)  <  yQ(t)  <  yo(T) 

and  h(x.,y.;  x  ;y  )  =0  elsewhere. 

1  J  1  o  o 

From  above  discussion,  we  can  find  a  solution  to  determine  the 
blurring  produced  by  linear  camera  motion.  We  will  again  describe 
the  image  g(xity^)  which  is  integrated  by  the  camera  sensor, 
as  follows 


g(x,y) 


T 

I  f[x  -  A(t) ,  y  -  B(t)  ]  dt 
n  o  o 


(2; 


Where  A(t)  and  B(t)  are  the  motion  in  x  and  y  directions 
respectively.  If  we  assume  that  camera  is  moved  0  degrees  off 
the  horizontal  axis  and  with  a  velocity  of  V  ,  then  the  displace¬ 


ment  distance  in  x  and  y  directions  at  an  instant  time  t  will  be 


N  % 


A(t)  =  V  t  CosQ 
B(t)  =  V  t  Sin0 


g(x,y)  =  /  f(x-V  tCosS,  y-V  t  Sin0)  dt 

0  c  C 


Taking  Fourier  transform  of  the  both  sides  of  Eq.  (30)  we  obtain 


T  00  00 

G(u,v)  =  /  /  /  f(x-V  t  Cos0,  y-V  t  Sin  ) 

0  c  c 


exp{  -  j2rr(ux+vy )}  dx  dy  dt 


Making  the  change  of  variables 


£  =  x-V  t  Cos0 

c 


n  =  y-V  t  Sin0 
c 


we  then  have 

T  00  CO 

G(u,v)  =  /  /  /  f(e,n)  exp{-j2f(ue+nv)}d£dri 

Q  _oo  _ oo 


exp{uV^tCos0  +  vV^t  Sin0}  dt 


=  F(u,v)  1  exp  {  -j2iT(uV  t  CosQ  +  vV  t  SinQ)}  dt  (33) 

0  c  c 

If  we  let  d  be  the  total  displacement  made  by  the  camera  in  the 
interval  (0,T),  then 


Crt’iTfriTC 


and  if  we  define  "f"  as 


f  =  u  Cos0  +  v  Sin6 


then,  Eq.  (33)  becomes 


G(u,v)  =  F(u,v)  /  exp{-j27rf  —  t}  dt 


as  we  can  see,  the  OTF  of  the  degradation  equals  to 


H  ,  ( u ,  v )  =  S  exp{-j2TTf^t}  dt 
0  1 


Performing  the  integration  in  Eq.  (37),  we  have 


Hd(u,v)  =  Sin  Trfd  e  jTTfd 


or 


Hd(u,v)  =  T  Sinc(fd)  e 
SinTrfd 


-jTTfd 


where  Sinc(fd)  = 


Trfd 


(35) 


(36) 


(37) 


(38) 


(39) 


OTF  of  Out-of-Focus  Lens  System 

An  aberration-free  (diffraction  limited)  imaging  system  converts 
a  diverging  spherical  wave  into  a  spherical  wave,  that  converges  toward 
a  point  in  the  image  plane.  An  aberration,  such  as  defocusing,  causes 
the  exit  wave  to  depart  from  its  ideal  spherical  shape.  In  Figure  3, 


0  is  an  axial  object  point  in  the  object  plane  and  I  is  the  geometrical 
image  of  that  point  in  the  image  plane.  The  energy  radiated  from  the 
point  object  0  passes  through  the  entrance  and  exit  pupils  and  makes 
an  angle  A  with  the  optical  axis  in  the  image  plane.  The  references 
s.  and  s.  are  centered  at  0  and  I,  respectively,  and  are  of  such 
radii  that  they  intersect  the  optical  axis  at  the  entrance  and 
exit  pupils,  respectively. 


Figure  3  Imaging  System  Without  Defect  of  Focus  (Ref.  8) 


Figure  4  Imaging  System  with  Defect  of  Focus  (Ref.  8) 

In  Figure  4,  s'  is  the  converging  wave  front  exiting  the  exit 
pupil  with  a  different  spherical  shape  from  its  ideal  shape  and  forms 
the  image  at  the  out-of-focus  point  P.  The  wave  front  converges 
to  the  in-focus  point  I.  We  can  specify  the  amount  of  defocusing 
by  the  out-of-focus  distance  Z  or  the  optical  distance  W  between 
S£  and  s'. 

The  pupil  function  which  describes  the  defocusing  effect  on 
the  ideal  spherical  wave  front  is  of  the  form  (Ref.  8) 


2  y 

'exp { jkw(x^  +  y  )} 


P(x,y)  = 


2  2 

if  x  +  y  <  1 


(40) 


where  k  =  2 tt/  X  and  W  is  the  optical  distance  between  the  emergent 

wave  front  and  the  ideal  spherical  wave  front  measured  along  the 

extreme  ray,  and  a  convenient  measure  of  the  defect  of  focus.  A 

is  the  wave  length.  The  coordinates  (x,y)  at  exit  pupil  are 

2  2 

normalized  so  that  x  +  y  =  1  at  the  edge  of  the  pupil. 

If  the  system  has  a  circular  aperture,  then  the  OTF  for 
such  a  system  is  a  function  of  only  one  spatial  frequency 
variable 

oo  oo 

/  /  P(x  +  y)  P*(x  -  y)  dx  dy 

_  CO  —  OO 

H(s)  =  - — -  (41) 

/  /  |P( x , y )  |  2  dx  dy 

_  OO  _ OO 

The  frequency  s  is  a  reduced  frequency  and  related  to 

2  2~ 

normalized  frequency  f  =  /  u  +  v  by 

s  -  (rk*)  f  (42> 

where  n  is  the  refractive  index  and  taken  to  be  1  for  all 
calculations  in  this  study  and  the  angle  A  is  defined  in 
Figure  4.  The  integrand  in  Equation  (41)  vanishes  for  s 
values  outside  the  range  0  £|  s|  £  2.  The  OTF  for  a  defocused 
system  can  be  found  by  inserting  the  pupil  function  of  Eq .  (40) 
into  Eq.  (41) 


Vs)  - 


—  Cos(ias)  fbJ.(a)  +  V  (-1)  Sm 
it  a  v.  1  , L . 

k=l 


4kb  T  (a)  T  (c 
2k  2k-l  2kn 


4  rl  ■,  r  /  ^k  Sin(2k+l)b 

-  -  sln  [j as]  J  (-i)  ^ 

k=0 


j  (a)  (a) 

X  J2k  2(k+l) 


where 


b  =  Cos  (s/2) 


and  J.  is  a  Bessel  function  of  first  kind  order  k.  These  series 
k 


are  slowly  converging.  The  OTF,  H_(s),  is  the  exact  or 
diffraction  OTF.  Normalization  in  Eq.  (40)  assures  that  Hj,(0)  = 


and  as  mentioned  before,  H^(s)  =  0  for  s  >  2.  Because  of  the 


symmetry  of  the  system,  Hg(s)  =  -Hg(s).  From  the  symmetry 


properties  of  the  Bessel  functions,  HP(s)  is  an  even  function  of 


the  defect  of  focus  distance  w. 


The  exact  OTF  obtained  in  Eq.  (43)  is  difficult  to  calculate 


but,  for  small  amounts  of  defocusing  and  small  angles  A,  this 


function  can  be  approximated  to  the  geometrical  OTF.  The 


geometrical  PSF  of  a  defocused  aberration  free  system  with  a 


circular  aperture  is  simply  a  circular  path  of  light  of  radius  R 


(Ref.  14),  that  is 


2  2 

where  R  =  /  x  +  y  .  The  fourier  transform  of  Eq.  (45)  will 
give  us  the  geometrical  OTF  of  the  defocused  system. 


Hd(f) 


J^(27trf  ) 


It  is  convenient  to  express  the  OTF  in  terms  of  the  optical  path 
length  W  and  the  normalized  frequency  s,  because,  the  OTF  can  then 
be  calculated  without  knowledge  of  the  optical  parameters  of  the 
system.  A  relation  between  W  and  the  out-of-focus  distance  Z  is 
required  to  calculate  the  OTF  for  a  specific  system.  From  the 
triangle  IKP  in  Figure  4,  we  can  write 

(KP)2  =  (KI)2  +  (IP)2  -  2( AI ) ( IP)  Cos  (tt-A)  (47 


This  can  be  written  as 


(r'+Z)^  =  (r'+w)2  +  Z2  -  2(r'+w)  Z  Cos  (tt-A) 


where  the  image  distance  01  =  r'.  Solving  for  W  gives 


W  =  -r’  -  Z  Cos  A  +  (r'2  +  2r'Z  +  Z2Cos  A)* 


This  equation  can  be  approximated  for  small  angles  and  small 


s  *. 
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W  =  Z  (1  -  Cos  A) 


(50) 


A  1 

By  using  the  trigonometric  small  angle  approximation  Sin  -j  =  —  Sin  A, 


we  find 


W  =  i  Z  Sin2A 


(51) 


We  can  again  approximate  the  argument  x  of  Bessel  function  in  Eq.  (46) 
by  using  Equations  (41)  and  (51) 


x  =  27irf 

=  Air  ws/  A  CosA 
=  (4tr/A)  ws 

x  =  a 


(52) 


Therefore , 

J,  (x) 

Hd(s)  .  2  —  (53) 

This  equation  should  be  normalized  so  that  H^(0)  =  1,  and  represents 
the  approximate  OTF  for  a  defocused  optical  system  with  circular 
aperture . 

Noise  Process 

Observations  of  an  image  field  are  subject  to  measurement  errors, 
irregularities  of  the  recording  media,  and  atmospheric  light 
fluctuations.  These  phenomena  give  rise  to  a  noise  in  the  image. 

Noise  in  images  can  be  divided  into  three  categories:  film  grain 


.  V 
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noise,  photoelectronic  noise,  and  thermal  noise. 


Film-grain  noise  is  produced  by  a  photographic  emulsion  during 
the  process  of  image  recording  and  reproduction.  Under  microscopic 
examination,  the  smooth  tones  of  a  photographic  image  assume  a 
random,  grainy  appearance.  Randomness  is  further  introduced  by  the 
variable  number  of  photons  which  expose  a  particular  film  grain  and 
the  varying  size  of  the  grains  themselves.  The  subjective 
appearance  of  these  factors  in  the  image  is  called  film-grain 
noise.  For  most  practical  purposes,  film-grain  noise  can  be 
modeled  a  Gaussian  process  (Ref.  4). 

Photoelectronic  noise  is  due  to  the  statistical  nature  of 
light  and  of  the  photoelectronic  conversion  process  inherent 
in  image  sensors.  A  model  of  this  noise  is  given  in  Ref.  5. 

Electronic  noise  is  due  to  random  thermal  motion  of  electrons 
in  resistive  circuit  elements  and  can  be  modeled  as  white  Gaussian 
noise  process;  i.e.,  it  has  a  flat  power  spectrum  (Ref.  5.). 

To  restore  a  picture  in  the  presence  of  noise,  in  addition  to  a 
knowledge  of  the  degrading  optical  transfer  function,  we  need  to  know 
(at  least  in  theory)  both  the  statistical  properties  of  the  noise  and 
how  it  is  correlated  with  the  original  image.  In  practice,  the  most 
common  assumptions  (like  in  the  electronic  communications)  about  the  noise 
are  that  it  is  white,  i.e.,  its  power  spectrum  density  is  constant,  and 
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it  is  uncorrelated  with  the  original  image  (signal).  Both  these 
assumptions  may  be  questioned.  The  concept  of  white  noise  is  a 
mathematical  abstraction,  however,  it  is  a  convenient  and 
reasonably  accurate  model,  if  the  power  spectrum  density  of  the 
original  image  falls  off  much  faster  than  the  power  spectrum  density 
of  the  noise.  As  far  as  the  noise  and  the  original  image  being 
uncorrelated  is  concerned,  there  are  examples  of  the  noise  which  is 
correlated  with  the  original  image,  e.g.,  in  the  case  of  film-grain 
noise  (Ref.  4,  5).  We  will  be  mainly  concerned  with  the  noise  which 
is  uncorrelated  with  the  original  image  and  additive  in  the  spatial 
frequency  domain.  Concerning  the  power  spectrum  of  the  noise,  we 
will  consider  both  cases  where  the  noise  is  white  and  where  it  is 
colored  in  the  spatial  frequency  domain.  The  latter  will  be 
modeled  as  a  gaussian  random  variable  with  zero  mean. 

Different  restoration  methods  may  require  different  amounts  of 
information  about  the  noise  process.  Constraint  image  restoration 
methods  (Ref.  7,  6)  require  that  only  the  variance  of  the  noise  be 
known.  On  the  other  hand,  the  Wiener  filter  and  power  spectrum 
equalization  filters  require  the  characterization  of  the  process  in 
terms  of  its  power  spectral  density  or  the  ratio  of  the  noise  power 
spectrum  density  to  the  original  image  power  spectrum  density  which 
we  refer  to  as  the  noise-to-signal  ratio  (NSR).  Since  we  will  use 
Wiener  filter  and  power  spectrum  equalization  filter  for 
restoration,  we  will  estimate  the  NSR  in  spatial  domain  and  in  the 


spatial  frequency  domain  by  assuming  that  the  noise  is  additive  and 
it  is  spatially  stationary  stochastic  random  process  with  zero  mean. 


Restoration  Methods 

As  we  mentioned  earlier,  image  restoration  is  an  estimation 
process  that  attempts  to  recover  an  original  image  from  its 
degradations  that  were  incurred  while  the  image  was  being  obtained. 
Restoration  techniques  are  methods  of  obtaining  a  PSF  (point  spread 
function)  for  the  restoration  filter,  so  that  the  output  of  the 
restoration  filter  is  a  good  estimate  of  the  original  image.  We  will 
assume  that  the  degraded  image  g(x,y)  and  the  original  image  f(x,y) 
obey  the  model  given  in  Eq.  (13).  That  is, 

OO  00 

g(x,y)  =  /  /  h(x-y,  v-8)  f (y,8)  dY  dB  +  n(x,y)  (54) 

— OO  — OO 


and  Fourier  transform  of  degraded  image  and  restored  image  are  given 
by  Equations  (15)  and  (16)  respectively. 


G(u,v)  =  Hd(u,v)  F(u,v)  +  N(u,v) 


F(u, v)  =  Hr(u,v)  G(u, v) 


(55) 

(56) 


Inverse  Filter. 

In  the  absence  of  noise.  Equation  (55)  reduces  to 


G(u , v )  =  H,(u,v)  F(u, v) 


(57) 


or  equivalently 


F(u , v) 

A 

since  we  want  F(u,v)  = 


=  G(u,v) 

Hd(u,v) 

F(u,v),  Eq.  (58)  implies  that 


(58) 


Hr(u,v) 


1 _ 

Hd(u, v) 


(59) 


and 


F(u,v) 


Hr(u,v)  G(u , v) 


G(u.v) 

Hd(u, v) 


(60) 


This  means  that  we  can  restore  the  degraded  image  g(x,y)  by 
multiplying  the  Fourier  transform  G(u,v)  of  the  degraded  image 
by  Hr(u,v)  and  then  inv  rse  Fourier  transforming. 

Some  problems  arise  when  one  wants  to  use  Eq.  (60)  in  practice. 
There  may  be  points  or  regions  in  the  Fourier  domain  where 
Hd(u,v)  =  0.  In  the  absence  of  noise,  as  we  can  see  from  Eq.  (57), 
the  Fourier  transform  G(u,v)  of  the  degraded  image  will  also  be  zero, 
leading  to  undetermined  ratios.  Therefore,  we  can  say  that  even  in  the 
absence  of  noise,  it  is,  in  general,  impossible  to  restore  f(x,y) 
exactly,  if  Hd(u,v)  has  zeros  in  the  Fourier  domain.  In  the  presence 
of  noise,  the  zeros  of  G(u,v)  and  Hd(u,v)  will  no  longer  coincide. 

When  noise  is  present,  we  have 


G(u , v )  =  H,(u,v)  F(u, v)  +  N(u , v) 


(61) 


applying  the  restoration  filter  gives 


F(u , v)  =  G(u,v)  Hr(u,v) 


F(u,v) 


N'(u,v) 
Hd(u, v) 


(62) 


If  Hd(u,v)  is  small  or  zero  in  any  region  of  the  Fourier  domain 
where  N(u,v)  is  not  correspondingly  smaller  or  where  N(u,v)  is 
non  zero,  then  the  contribution  of  the  noise  will  be  amplified,  thus 
the  term  \( u , v  ) , F j( u , v )  will  have  much  larger  magnitude  than  F(u,v) 
and  the  inverse  Fourier  transform  of  G(u , v )H  ( u , v)  will  be  strongly 
influenced  by  these  large  terms.  Instead,  it  will  contain  many 
noiselike  variations  and  will  not  be  meaningful  restoration  of 
f(x,y). 

Since  the  right-hand  side  of  Eq.  (60)  cannot  normally  be 
evaluated  numerically  due  to  the  zeros  of  H^Cu.v),  the  most  we  can  do 
is  to  restore  only  the  frequencies  in  the  Fourier  domain  where  the 
noise-to-signal  ratio  is  low.  Such  an  alteration  in  the  restoration 
process  means  that  the  exact  inverse  filter  is  not  being  performed; 
instead,  an  approximation  is  being  made.  This  approximation  can  be 
improved  to  some  degree  by  examining  the  output  of  the  restoration 
filter.  Let  us  denote  the  approximated  inverse  filter  transfer 
function  by  A(u,v) .  We  can  say  that 

Hr(u,v)  =  A(u, v)  (63) 

If  H,(u,v)  does  not  have  zeros  and  noise  is  not  present,  then  the 


restoration  filter  optical  transfer  function  will  be 

Hr(u,v)  =  A(u,v)  =  H  ^  ys)  (64) 

Wiener  Filter 

The  inverse  filter  performs  poorly  in  the  presence  of  noise  since 
the  noise  process  is  not  taken  into  account  in  the  filter  design. 
Inverse  filtering  also  requires  tedious  manipulation  of  the  inverse 
of  the  degrading  transfer  function.  It  would  be  preferable,  if  the 
filter  itself  controlled  the  noise  amplification  automatically.  One 
way  to  design  such  a  filter  is  to  find  a  restoration  filter  transfer 
function  that  minimizes  some  measure  of  the  difference  between 
estimated  image  f(x,y)  and  the  original  image  f(x,y).  We  will  design 
a  restoration  filter  that  minimizes  (in  a  statistical  sense)  the  mean- 
squared  error  between  the  original  image  f(x,y)  and  the  estimated 
image  f(x,y).  This  filter  is  called  least  squares  filter  or  Wiener 
filter . 

Let  the  original  image,  the  degraded  image  and  the  noise  be 
represented  by  the  stochastic  random  processes  f(r),  g(r),  and  n(r) 
respectively.  Where  r  is  a  position  in  the  xy-plane.  Then  from  the 
superposition  integral  in  Eq .  (9),  we  can  write 


"Vr '  \\\  v v 
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where  dr'  represents  the  area  element  dx’dv'  in  the  xy-plane  and  h^(r) 
is  the  point  spread  function  of  the  degradation.  We  need  to  solve  the 
integral  Equation  (65)  for  the  original  image  process  f(r).  However, 
with  n(r)  unknown,  the  Eq .  (65)  cannot  be  solved  directly.  The  most 
we  can  do  is  to  produce  an  estimate  f(r)  of  the  solution.  We  assume 
that  f(r)  and  n(r)  can  be  viewed  as  spatial  stochastic  processes,  and 
assume  that  noise  process  n(r)  is  known.  Then,  the  best  estimate  is 
to  maximize  the  posterior  probability  density  of  f(r),  given  g(r),  as 
defined  Bayes's  rule  (Ref.  21).  If  we  further  assume  that  the 
estimate  f(r)  is  a  linear  function  of  the  gray  levels  in  the  degraded 
image  process  g(r),  then  the  linear  estimate  f(r)  can  be  expressed  as 


A 

f(r)  =  :  hr(r-r’)g(r')  dr'  (66) 

Our  aim  is  to  find  a  function  hr(r),  which  is  the  point  spread 
function  of  the  restoration  filter,  such  that  mean-square  error 


£  2  =  E  {  [  f  ( r )  -  f(r)]2} 

=  E  { [f (r)  -  f  hr(r-r')  g(r’)  dr’  ]2} 


(67) 


is  minimized. 

A  function  which  satisfies 


E{g(s)[f(r)  -Jhr(r-r')  g(r')  dr']}  =  0 


(68) 


for  all  positions  r  and  s  in  the  xy-plane  will  also  minimize 
Eq .  (67).  From  Eq.  (68),  we  can  write 
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By  using  the  definitions  of  the  autocorrelation  and  cross-correlation 
of  a  random  process  (Ref.  9),  we  can  re-write  Eq.  (69)  as 


m  ww  m  wrwr*:  w;  wvw.w*.  "jcw  v  v  pji  v  v  ■>  ■;  vj.'  wum  w»  v*  v»v»v^  v*  ^  w.* 


or 


v  •■• 


h  ( x , y )  *  R  (Ax,  Ay)  =  R  f(Ax,  Ay) 
re  et 


(75) 


If  we  take  the  Fourier  transform  of  both  sides  of  Eq .  (75)  we  have 


H  ( u , v )  S  ( u  ,  v )  =  S  f(u,v) 


(76) 


or 


H  (u.v) 

l 


S  -  ( u ,  v  ) 

_ai _ 

Sg(u.v) 


(77) 


where  S  ,(u,v)  and  S  (u.v)  are  the  Fourier  transforms  of  R  f(u,v) 


gf 


gf 


and  R  (u,v),  respectively,  R  ^(Ax,  Ay)  is  the  spatial  cross- 
g  g 

correlation  of  the  degraded  image  g(x,y)  and  original  image  f(x,y) 
and  it  is  given  by 


R  (Ax,  Ay)  =  E(g(x,v)  f(x+Ax,  y+Av)} 

gt 

R  f .Ax,  Ay)  =  E{[hd(x,y)  *  f(x,y)  +  n(x,y)] 

f(x+Ax,  y+Ay)} 


(78) 


(79) 


Since  we  assume  that  image  process  and  noise  process  are  uncorrelated 
and  that  the  noise  process  has  zero-mean,  we  have 

E(n(x,y)  f(x+Ax,  y+Ay)}  =  E{n(x,y)}  E{f(x+Ax,  y+Ay)}=  0  (80) 


and  Equation  (79)  becomes 


R  j (Ax,  Ay)  =  hd(x,y)  *  Rf(Ax,  Ay) 


(81) 
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Taking  the  Fourier  transform  of  the  both  sides  of  Eq.  (81)  gives  us 

S  ,(u,v)  =  H*(u, v)  Sf(u,v)  (82) 

g-  d  t 

where  Sr(u,v)  is  the  power  spectrum  of  the  original  image.  It 

can  easily  be  seen  from  Figure  2  that  S  (u,v)  is  given  by 

8 

Sg(u, v)  =  |Hd(u,v)|  2  Sf(u,v)  +  Sn(u,v)  (83) 

where  sn(u,v)  is  the  power  spectrum  of  the  noise  process.  If  we 
insert  Equations  (82)  and  (83)  into  Equation  (77),  we  obtain 


H*(U,V)  Sr (U , V) 

H(u,v)  =  - - - - 

|Hd(u,v)f  Sf(u,v)  +  Sn(u,v) 


Eq.  (84)  can  be  written  in  more  convenient  form  as 


(84) 


Hr(u,v) 


H*(u,v) 
|Hd(u,v)|2  + 


Sn(u, v) 

Sf(u, v) 


(85) 


This  result  gives  the  transfer  function  of  Wiener  filter  that  we 
will  use  as  the  restoration  filter  transfer  function. 


Power  Spectrum  Equalization  Filter. 

We  modeled  the  degraded  image  as  the  output  of  a  linear  system 
which  has  the  degrading  transfer  optical  function  Hd(u,v).  The  input 
to  this  system  is  taken  to  be  original  image.  From  linear  systems  theory, 
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we  can  establish  a  relationship  between  the  original  image  power 
spectrum  and  the  degraded  image  power  spectrum.  From  Figure  2 
the  power  spectrum  of  the  degraded  image  can  be  written  as 

S  (u,v)  =  |  Hd(u , v)|  2  Sf(u,v)  +  Sn(u,v)  (86) 

Power  spectrum  equalization  criterion  states  that  the  power  spectrum 
of  the  estimated  image  be  equal  to  the  power  spectrum  of  the  original 
image.  That  is, 


Sf(u,v)  =  S*(u,v) 


(87) 


where  S^(u,v)  denotes  estimated  image 
possible  to  write  S^(u,v)  in  terms  of 
function  of  the  restoration  filter 


power  spectrum.  It  is 


S  (u,v)  and  the  transfer 

g 


Sj(u,v) 


|  Hr(u,v)  ]2  Sg(u, v) 

|Hr(u,v)|2  [|Hd(u,v)  2  Sf(u,v)  +  Sn(u,v)|]  (88) 


by  substituting  Equation  (87)  into  Equation  (88)  and  solving  for  Hr(u,v) 
we  obtain  the  transfer  function  of  the  restoration  filter  that  performs 
the  power  spectral  equalization  as  follows 


Hr(u,v) 


Sn(u, v) 
+  Sf(u,v) 


1 

2 


(88.1) 
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III.  Parameter  Estimation  Techniques 


Estimation  of  Degrading  Transfer  Function 

We  have  seen  several  methods  for  restoring  degraded  images. 

These  methods  require,  as  a  description  of  the  degradation,  the 
degrading  point  spread  function  (PSF)  or  its  Fourier  transform,  the 
degrading  optical  transfer  function  (OTF).  However,  in  many  practical 
situations,  a  prior  knowledge  of  the  degrading  PSF  or  OTF  is  not 
available.  Therefore,  in  order  to  apply  one  of  these  restoration 
methods,  one  must  first  estimate  the  degrading  OTF.  Since  we  assume 
that  we  do  not  have  a  prior  information  about  the  degradation  and  that 
the  original  undegraded  image  is  unavailable,  we  will  consider  the 
estimation  of  the  degrading  OTF  from  the  degraded  image  itself;  i.e., 
so  called  "blind"  restoration. 

Method  of  Frequency  Domain  Inspection. 

The  degraded  image  in  the  absence  of  noise  can  be  expressed  as  in 
Eq.  (57)  in  the  frequency  domain.  It  is  given  by  point-by-point 
product  of  the  degrading  OTF  H^Cu.v)  and  the  Fourier  transform  of  the 
original  image.  Therefore,  on  the  lines  where  H^(u,v)  is  zero,  G(u,v) 
will  also  be  zero.  If  magnitude  of  G(u,v)  is  displayed,  these  lines  of 
zero  values  can  be  determined,  given  that  H^Cu.v)  is  not  zero  in  the 
neighboring  areas  (Ref.  10).  In  the  presence  of  noise,  the  displayed 


values  may  not  be  exactly  zero  at  these  points,  but  the  rough 
minima  corresponding  to  ideal  zeros  can  still  be  identified,  if  the 
noise  is  sufficiently  small  or  white. 

The  identification  of  the  zeros  in  the  OTF  may  not  be  sufficient 
to  determine  the  complete  function.  However,  many  of  the  common  causes 
of  the  degradation  such  as  the  blurs  due  to  a  relative  motion  between 
camera  and  object  or  the  out-of-focus  camera  have  distinct  zero  patterns  in 
their  OTF's.  We  will  ignore  the  unlikely  possibility  that  the  zeros  of 
the  Fourier  transform  of  the  original  image  have  the  same  zero  pattern, 
we  will  also  assume  that  the  degradation  is  due  to  either  motion  blur 
or  out-of-focus  camera.  If  the  restored  image  looks  reasonable,  the 
confidence  in  the  correctness  of  these  assumptions  increases.  In  the 
case  t..at  some  prior  knowledge  of  the  OTF  is  available,  this  knowledge 
can  be  used  in  combination  with  the  above  method.  For  example,  if  the 
shape  of  the  OTF  is  known,  then  only  the  frequency  scale  factor 
needs  to  be  determined.  In  fact,  in  some  degraded  images,  it  is 
possible  to  determine  the  degradation  form  interactively.  Once  correct 
type  of  OTF  is  determined,  the  accuracy  of  its  scale  factor,  in  some 
cases,  can  be  improved  by  averaging  J G( u , v ) [  or  |G(u,v)|^  along  the 
lines  parallel  to  the  expected  zeros  and  examining  the  resulting 
function . 

In  using  this  method,  it  is  desirable  to  perform  a  simple  linear 
spatial  filtering  operation  such  as  differentiation  on  the  degraded 
image  before  Fourier  transforming  it,  in  order  to  suppress  the  low 


frequencies  and  to  enhance  the  high  frequencies.  In  most  images, 
the  amplitudes  of  the  low  frequency  components  are  much  greater  than 
the  amplitudes  of  the  high  frequency  components.  The  above  filtering 
operation  may  flatten  the  frequency  spectrum  of  the  degraded  image 
and  makes  the  identification  of  the  zeros  easier. 

One  such  filter  which  meets  the  above  objective  is  Laplacian 
operator  7^  =3^/3x"  +  9  /3y“.  Since  the  Laplacian  operator  is 
linear  and  shift-invariant,  applying  it  to  the  picture  g(x,v)  is 
equivalent  to  convolving  the  degraded  image  g(x,y)  with  the  PSF  h,(x,y) 
defined  by  (Ref.  11). 

0  1  0 

hr(x,y)  =  1  -4  1  (89) 

0  10 

Visualization  of  the  zeros  in  the  spectrum  can  also  be  aided  by 
taking  the  logarithm  of  the  degraded  image  power  spectrum.  This 
operation  may  enhance  the  amplitude  of  dips  due  to  the  degrading 
transfer  function.  If  the  zeros  are  equally  spaced,  they  produce  a 
series  of  periodic  spikes  in  the  log  power  spectrum.  The  power 
spectrum  of  the  log  power  spectrum,  called  the  cepstrum,  can  be  used  for 
determining  the  exact  spacing  of  the  spikes  and  consequently  the  zeros 
of  the  degrading  transfer  function. 

After  visualization  of  the  zeros,  we  need  to  recognize  the  zero 
pattern  of  the  degrading  OTF.  In  the  case  of  linear  camera  motion, 
the  zeros  of  the  OTF  requires  a  linear  pattern.  This  pattern  is  the 


been  studied  by  many  researchers  (Ref.  12,  13,  14,  15),  is  as 
follows:  the  degraded  image  is  divided  into  M  regions  (subimages) 
which  are  all  identical  in  size  and  they  may  overlap.  If  we  assume 


r.*. 


-.*' 


that  the  extent  of  the  degrading  point  spread  function  h^(x,y)  is 
small  compared  to  that  of  each  region  of  the  degraded  image,  then 
approximately  we  have 

00  CO 

gi(x,y)  =  /  /  hd(x-y,  y-3)  fi(Y,S)  dY  d6  +  ni(x,y) 


(92) 


where  the  index  i  denotes  the  i'th  subimage.  Since  we  can  pre-process 
(low-pass  filter)  the  degraded  image  to  reduce  the  noise  effects,  we 
ignore  the  noise  term  in  the  right-hand  side  of  Eq.  (92).  By  taking 
the  Fourier  transform  of  the  both  sides  of  Eq.  (92),  we  obtain 

Gi(u,v)  =  Hd(u,v)  Fd(u,v)  ,  i  =  1,  2,  3 . M  (93) 

The  functions  Gd(u,v),  Fd(u,v),  Hd  (u,v)  are  complex  in  the  Fourier 
domain.  Therefore,  we  may  write  Eq.  (93)  as 


|Gi(u, v) | e 


j0j(u,v) 


=  ! Hd ( u , v ) |  ‘  |Fd(u,v) |  •  e 


j[0H(u,v)  +  0fi(u,v) ] 

(94) 


By  taking  the  product  average  of  the  magnitudes  over  M  regions,  we  will 
have 


M 

M 

n  I g  (u, v) 1  = 
i=l 

n  1 

i=l 

I  F.(u,v)| 

|Hd(u,v) | 

M 


it  is  possible  to  estimate  the  magnitude  of  the  degrading  OTF  by 
(Ref.  13) 


(95) 
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have  information  about  the  phase  of  the  estimated  OTF.  In  most  cases, 
it  might  be  reasonable  to  assume  that  the  phase  of 

M  1/M 

II  |  F.(u,v)|  is  approximately  constant. 

i=l  1 

Another  way  of  estimating  the  OTF  by  segmentation  is  to  treat  the 
image  as  a  (spatially  stationary  stochastic  process  (Ref.  14).  From 
the  stochastic  linear  system  relationship,  we  can  write  the  power 
spectrum  of  the  degraded  image  in  the  i'th  region  as 


S  •  ( u ,  v )  =  Sf.(u,v)  |  H,(u,v)|  1  +  S  .  ( u ,  v )  (97) 

gi  tx  1  a  1  ni 

After  taking  an  average  over  the  power  spectrum  of  the  subimages  (the 

square  of  the  magnitude  of  the  Fourier  transform  of  the  subimages),  we 

have  an  estimate  of  the  power  spectrum  of  subimages.  Such  an  average 

will  retain  much  of  the  flavor  of  |Hj(u,v)|,  since  |H^(u,v)|  is  included 

in  each  one  of  the  averaged  subimages,  while  each  F^(u,v)  and  N^(u,v) 

are  varying  from  subimage  to  subimage,  thus,  introducing  a  much  more 

subtle  contribution  to  the  average  (Ref.  14).  We  may  now  search  for 


the  zeros  of  H^(u,v)  more  successfully  than  we  did  in  the  frequency 
domain  inspection  method.  Since  most  degraded  images  have  reduced 
amount  of  high  frequency  power,  the  search  for  the  zeros  of  |H^(u,v)  | 
should  be  carried  out  along  the  dc  axes  of  estimated  S  (u,v).  The 

O 

power  cepstrum  which  is  defined  as  the  Fourier  transform  of  the 
logarithm  of  the  power  spectra  Se(u,v)  may  enable  us  to  identify 

O 

exact  spacing  of  the  zeros  of  H^(u,v).  The  power  cepstrum  of  Eq . 

(97)  is  equal  to 

P(u,v)  =  F{log  Sf  ( u ,  v )}  +  2  F{log  Hj(u,v)  | }  (98) 

Note  that  we  have  neglected  the  noise  term.  Since  the  right  hand 
side  of  Eq.  (98)  results  from  an  averaging  process,  F{log  Sr(u,v)} 
retains  little  of  the  flavor  of  the  original  image  while  2  F{log|H^(u,v 
remains  strongly  characteristic  of  the  degrading  OTF.  The  rings  of 
spikes,  as  previously  mentioned,  results  from  focus  blur  and  the 
radius  of  these  rings  will  be  the  parameter  that  we  need  to  know 
to  estimate  the  degrading  OTF.  In  the  case  of  motion  blur,  these 
spikes  will  be  periodic  in  the  direction  of  the  motion  and  their 
pattern  is  given  by  Eq.  (90). 

Another  approach  which  is  similar  to  the  method  which  is  defined 
in  Ref.  12  may  be  used  for  estimating  the  degrading  OTF  by  segmentation 
method.  By  taking  the  complex  logarithm  of  the  both  sides  of  Eq .  (93), 
we  have 


CL0G[G. (u, v) ]  =  CLOG [ F . ( u , v ) ]  +  CLOG  [H,(u,v)] 


(99) 


Averaging  over  M  subimages, 


M  M  M 

^  l  CL0G[G.(u,v)]  =  l  CLOG  [ F. (u , v) ]  +  ^  [  CLOG[  Hd(  u ,  v )  ]  (100) 

i=l  ‘  i=l  i=l 


We  assumed  that  H^(u,v)  is  included  in  each  subimage.  Therefore 
we  can  write  Eq.  (100)  as 


M  M 

CL0G[Hd(u,v)]  =  i  l  CL0G[G.(u,v)]  -  ±  l  CL0G[ F. (u, v) ]  (101) 

i=l  ‘  i=l 


Since  the  right  hand  side  of  Eq .  (100)  results  from  an  averaging  proces 

and  the  term  CL0G[H  ,(u, v) ]  is  included  by  each  subimage,  the  averaged 
1  M 

term  —  V  CL0G[ F. (u, v) ]  retains  little  of  the  flavor  of  CL0G[ F . ( u , v ) ] . 

•1.^1  L 

1  l 

If  the  term  —  )  CL0G[F. (u,v) ]  was  available,  one  could  estimate  the 
i=l  1 

complex  logarithm  of  the  degrading  0TF  by  Eq.  (101).  However,  we  must 
estimate  the  degrading  0TF  without  a  prior  knowledge  of  the  original 
image  and  the  degradation.  If  the  intensity  distributions  in  the 


M  regions  of  the  original  image  are  sufficiently  different  from  one 

another,  and  if  M  is  large  enough  (M  >100),  then  it  may  be  assumed  that 

1  ? 

the  term  —  l  CL0G[F . ( u, v) ]  converges  to  zero.  Therefore,  we  have 
M  i=l  1 


l 

CL0G[Hd(u,v)]  =  ^  [  CL0G[Gi(u,v)]  (102) 


Hd(u,v)  is  then  obtained  by  taking  the  anti logarithm  i.e., 


Hd(u,v)  =  exp  l  CLOG[G.(u,v)]^  (103) 

We  can  estimate  the  magnitude  of  the  degrading  OTF  by  one  of 
the  above  methods.  We  still  need  to  estimate  the  phase  of  the 
degrading  OTF.  A  method  for  this  estimation  is  given  in  Ref.  15  and 
will  be  summarized  here.  After  segmenting  the  degraded  image  and 
taking  the  Fourier  transform  of  each  segment,  we  have  (Ref.  Eq.  93) 

0Gi(u,v)  =  9 j_j ( u ,  v )  +  9fi(u,v)  (104) 

By  averaging  over  M  regions, 

M  M 

i  I  QG.(u,v)  =  9h(u,v)  +  ^  l  9f.(u,v)  (105) 

i=l  '  i=l 

Our  hope  was  that  the  second  term  in  the  right-hand  side  of  Eq .  (105) 
would  converge  to  a  constant  value.  However,  we  can  not  make  sure  that 
this  term  converges  to  anything  meaningful.  Another  method  for 
estimating  the  phase  (Ref.  15)  is  to  form  the  product  of  subimages 

G^(u,v)  G0:‘(u+Au,  v+Av)  =  H^(u,v)  H^’:‘  (u+Au,  v+Av)  F  (u,v) 

F.'::'(u+Au,  v+Av)  (10b) 

i 


and  taking  average  over  M  regions,  we  obtain 


'jr*y,*ymy  »j  ■?»>  vv?.^1 


v. 


1 


M 


—  y  Gi(u,v)  G^u+Au,  v+Av)  =  Hd(  u ,  v)Hd"(u+Au  ,  v+Av) 
i=l 


M 


4  I  F  (u,v)F  “(u-tAu,  v-tAv) 
i=l 

M 


(107) 


In  Equation  (107),  we  can  calculate  -^r  T  G.(u,v)  G.*(u  +Au,  v  +Av) 

n  i=l  1  1 


from  the  degraded  image  by  dividing  the  degraded  image  into  M  regions 
(subimages),  and  Fourier  tansforming  each  subimage  to  form  the  auto¬ 


correlation  product  G^(u+v)  G^Ku-fA  u,  v+Av)  (where  Au  and  Av  are 


constants  for  M  subimages,  and  then  averaging  the  resulting  products 


by  M.  If  a  prototype  image  which  is  similar  to  the  original  image  is 

1  ? 

available,  then  we  can  obtain  a  rough  estimate  of  77  l  F.(u,v) 


M  .  “  i- 
i=l 


Fd*(u+Au, v+Av)  from  the  prototype  image.  An  estimation  of  the  product 
Hd(u,v)  Hd*(u+Au,  v+Av)  can  then  be  made  by  Eq.  (107).  For  a 
lossless  degration,  which  means  the  volume  under  f(x,y)  is  equal 
to  the  volume  under  g(x,y),  we  have  (0,0)  =  1.  Assuming  that 
the  products  can  be  estimated  and  given  ( 0 , 0 )  =  1 ,  we  can  obtain 


a  recursive  relationship  where  H^^u+Au,  v+Av)  is  expressed  in  terms 


of  Hd(u,v) 


Hd‘”'(u+Au,  v+Av) 


~  i  G . (u , v )  G . *( u+Au ,  v+Av) 
i=l 


Hd(u,v)  yj  T  F^ ( u ,  v)  F^fu+Au,  v+Av) 
i=l 


M 


(108) 


However,  the  stability  of  this  expression  is  quite  poor.  For  the 
phase  of  0TF,  we  have 
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r. 


M 


A 


3 


-  A 


+  0p^(u,v)  -8p^(u+Au,  v-tA v)  (109) 


If  Eq.  (109)  is  averaged  in  some  sense  over  i,  and  if  the  average 

value  of  [0c.(u,v)  -  0..,.(u4Au,  v+Av)]  can  be  estimated,  then  noting 
II  r  1 

that  0  (0,0)  =  0,  we  may  estimate  0^(u,v)  as  (Ref.  15). 

0H(u,v)  =  ©H(u,v)  -  [0  (u,v)  -  6Gi(u4Au,  v+Av)]  av 

-  [6  p^(u,v)  -0p_^(u+Au,  v+Av)]  av  (110) 

Line  Estimation  of  the  OTF. 

Some  degrading  optical  transfer  functions  such  as  motion  blur 
in  one  direction  can  be  described  by  a  single  line  which  may  be 
estimated  from  the  Fourier  transform  of  the  degraded  image  itself. 

To  support  this  statement,  consider  the  degrading  OTF  of  the  motion 
blur  given  by  Eq.  (39). 

Hj(u,v)  =  T  Sinc(fd)  e 

where  f  =  u  Cos0  +  V  Sin0.  If  we  assume  that  the  relative  motion 
between  camera  and  object  is  in  the  horizontal  direction  (0=0),  the 
Equation  (39)  becomes 


Hd(u, v) 


(111) 


=  T  Sinc(ud)  e  '^7Tud 

or  for  the  motion  in  vertical  direction  (  0  =  -^) 

Hd(u,v)  =  T  Sinc(vd)  e  jTTvd  (112) 

Since  the  right-hand  side  of  the  Eq.  (Ill)  does  not  depend  on  the 
vertical  spatial  frequency  variable  v,  we  can  say  that 

Hd(u,v)  =  Hd(u)  =  T  Sine  (ud)  e  '*7Tud  (113) 

Recall  that  the  degraded  image  can  be  manipulated  by  a  digital  computer 
in  N  by  N  matrix  form  where  the  rows  and  columns  of  this  matrix 
represent  the  spatial  frequency  variable  v  and  u,  respectively.  In  the 
absence  of  noise,  the  Fourier  transform  of  the  degraded  image  is  equal 
to  the  point-by-point  multiplication  of  the  Fourier  transform  of  the 
original  image  and  the  degrading  OTF.  Therefore,  the  degrading  OTF 
and  the  original  image  matrices  has  to  be  identical  in  size.  We  can 
modify,  Eq.  (113)  to  write  it  in  matrix  form 

Hd(k,£)  =  T  Sinc(ud)  e~jTTud  (114) 

where  k,£  denotes  k'th  column  and  & ' th  row  of  the  degrading  OTF.  And 
if  the  degraded  image  is  represented  by  an  NxN  matrix,  k  and  l  take 
values  from  1  to  N.  If  the  function  given  by  Eq .  (114)  is  centered 
at  N/2,  then  we  can  say 


Equations  (113)  and  (114)  show  that  the  degrading  OTF,  in  this 
particular  case,  consists  of  identical  rows  (lines).  Therefore, 
it  will  be  sufficient  to  obtain  a  row  of  the  degrading  OTF. 

In  order  to  estimate  a  row  of  the  degrading  OTF,  we  will 
follow  a  procedure  similar  to  the  segmentation  method.  By  taking 
the  complex  logarithm  of  the  Fourier  transform  of  the  degraded  image, 
we  obtain 

CL0G[G(k,4)]  =  CL0G[Hd(k3l  )]  +  CL0G[F(k,4)]  (116) 

where  k,4  =1,  2,  3,  ....  N.  Averaging  over  N  rows,  we  have 

N  N 

~  V  CLOG[G(k , 1 ) ]  =  h  l  CLOG[H  ,  (k,  4)  ] 

4=1  ‘  4  =  1 

N 

+  4  T  CLOG[ F(k,4  ) ]  (117) 

4  =  1 

We  can  write  from  Eq .  (117) 

N  N 

|  I  CL0G[G(k,4  )  ]  =  CLOG[H  (k)  ]  +  4  I  CL0G[F(k,4)]  (118) 

‘  4=1  a  N  4=1 

In  order  to  obtain  an  estimate  of  the  first  term,  CLOG[Hd(k)], 

in  the  right-hand  side  of  Eq.  (118),  we  need  to  know  the  second 
1  N 

term,  tj-  T  CLOGf  F(k  X  )  ] .  Since  right  hand-side  of  Eq .  (118)  results 
4=1 

from  an  averaging  process,  and  the  term  CLOG[Hd(k)]  is  constant  while  the 
term  CL0G[F(k,4)]  is  changing  randomly  from  pixel  to  pixel  on  the  same 

1  N 

line  (4  th  row),  the  resultant  average  term  tt  V  CL0G[G(k,4)]  retains 

N  4  =i 


1  1V 

little  of  the  flavor  of  CLOG[F(k,£  ) ] .  If  the  term  ±  l  CLOG[F(k,£)] 


can  be  neglected,  then  we  have 

1  N 

CLOG [  H  ,  ( k )  =  ±  l  CLOG[G(k, £) ] 

d  N  £=1 


H.(k)  =  exp  ~  l  CLOG[G(k,£  )  ] 

a  N  £  =1 


(119) 


(120) 


If  the  light  distribution  in  the  degraded  image  varies  sufficiently 
fast,  we  can  expect  that  the  second  term  in  the  right-hand  side  of 
Eq.  (118)  can  be  approximated  by  a  constant.  In  this  case,  we 


obtain 


1  * 

H.(k)  =  exp  £  l  CLOG  [G(k,£)]  -  c 


(121) 


where 


C  I  CLOG  [  F(  k  ,£  )  ] 


(122) 


A  prototype  image  which  is  similar  to  degraded  image  can  also  be  used 
1  N 

for  estimating  tt  £  CL0G[F(k,£ ) ] .  Let  P(k,l)  be  the  Fourier 
L=1 

transform  of  the  prototype  image,  then  we  may  assume  that 


-i  N  ,  N 

77  I  CL0G[F(k,£)]  %  [  CL0G[P(k,£)] 

£=1  N  =1 


(123) 


Since  the  degrading  OTF  ought  to  be  identical  in  size  to  the  degraded 
image,  we  need  to  obtain  a  degrading  OTF  which  is  in  the  form  of  NxN 
matrix  from  the  estimated  row,  Hj(k).  This  can  be  done  by  duplicating 
H^(k)  N  times. 


48 


We  have  discussed  another  approach  called  "Line  estimation 
of  the  OTF"  to  estimating  the  degrading  OTF  of  the  motion  blur  in  one 
direction.  It  seems  that  it  is  possible  to  estimate  the  OTF  of  such 
a  degradation  from  the  degraded  image  itself  by  this  method.  The 
only  problem  in  this  approach  seems  to  be  the  estimation  of  the  term 

1  N 

—  £  CLOG[F(k.Jl )  ] .  We  consider  two  cases  to  overcome  this  problem. 

Z=1 

We  first  assumed  that  this  term  would  converge  to  zero  or  a  constant, 
and  secondly  we  estimated  it  from  a  prototype  image  which  is  assumed 
to  be  similar  to  that  of  original  image.  However,  we  can  say  that  both 
cases  may  contribute  an  error  to  estimated  degrading  OTF  due  to  the 
1  N 

estimation  of  -n-  £  CLOG[F(k,&) ] .  In  order  to  avoid  the  estimation 
i=l 

1  N 

of  the  term—  £  CLOG[F(k,£) ] ,  we  will  employ  another  approach  which 
i=l 

is  similar  to  the  frequency  domain  inspection  method.  Since  the  first 
term  in  the  right-hand  side  of  Eq.  (118)  converges  to  the  complex  logarithm 
of  the  degrading  OTF,  it  may  be  possible  to  recognize  the  zero  pattern 
of  the  degrading  OTF  correctly.  The  only  parameter  we  need  to  know 
then  is  the  total  displacement  distance  d  in  Eq .  (114).  If  the  first 
zero  of  the  pattern  is  "A"  pixels  (columns)  away  from  the  origin,  then 
the  total  displacement  distance  d  will  be 

d  =  j  (124) 

After  determining  the  parameter  d  of  the  Eq .  (124),  we  can  easily 


calculate  Eq.  (114)  for  a  given  exposure  time  T. 


Logarithmic  Estimate. 

In  cases  where  the  previous  methods  are  not  successful  in 
estimating  the  degrading  OTF,  we  can  obtain  an  approximation  by  the 
following  method.  In  the  absence  of  noise,  the  complex  logarithm 
of  the  Fourier  transform  of  the  degraded  image  is  given  by  Eq.  (116). 
For  an  NxN  degraded  image,  we  divide  both  sides  of  Equations  (116)  by 
N  to  obtain 

^  CL0G[G(k.O]  =  |  CL0G[Hd(k,5.)]  +  ^  CLOG[F(k,£)]  (125) 

Since  Hd(k.£)  is  spatially  invariant,  the  first  term  in  the  right-hand 
side  of  Eq.  (125)  will  be  the  logarithm  of  the  degrading  OTF  scaled  by 
1/N,  while  the  second  term,  -jjj-  CL0G[F(k£  )],  converges  to  a  small 
number.  Therefore,  we  may  say  that  the  right-hand  side  of  Eq.  (125) 
resembles  the  degrading  OTF  rather  than  the  original  image;  i.e.,  a 
rough  approximation  of  the  degrading  OTF  is 

Hd(k,£)  =  |  CLOG[G(k,£)]  (126) 


Estimation  of  Noise-to-Signal  Ratio. 

The  least-squares  and  the  power  spectrum  equalization  methods 
require  the  noise-to-signal  power  density  ratio  in  the  degraded  image 
as  a  parameter.  We  can  also  say  that  these  methods  respond  only  to  the 


ratio  of  the  noise-to-signal  power  density.  Since  we  do  not  have  a 
prior  knowledge  about  the  original  image  and  the  noise,  we  must  estimate 
the  NSR  (noise-to-signal  power  ratio)  from  the  degraded  image  itself. 


It  is  also  possible  to  define  the  NSR  as  the  variance  of  the  noise 
process  (Ref.  18),  i.e. 

NSR  =2 a  2  (130) 

n 

We  can  apply  one  of  above  equations  to  estimate  the  NSR;  however, 
we  first  need  to  estimate  the  variance  of  the  noise  process  and 
the  original  image.  Consider  a  relatively  noisy  region  in  the  degraded 
image.  This  region  consists  of  the  blurred  image  and  additive  noise. 


Since  the  region  is  relatively  noisy;  i.e.,  it  includes  more  noise 
components  than  the  blurred  image  components,  the  variance  of  this 
region  may  resemble  the  variance  of  the  noise  rather  than  the 
variance  of  the  original  image.  Therefore,  it  is  reasonable  to  assume 
that  the  variance  of  a  relatively  noisy  region  in  the  degraded  image 
be  approximately  equal  to  the  variance  of  the  noise  process.  If  we 
pre-process  (high  pass  filter)  the  degrading  to  enhance  the  high 
frequency  components  than  the  original  image  does,  and  if  the  variance 
is  calculated  over  several  noisy  regions  in  the  degraded  image  and  the 
resulting  variances  are  averaged  to  obtain  the  noise  variance,  then 
the  confidence  in  above  assumption  increases.  The  Laplacian  operator 
may  satisfy  the  objective  of  above  pre-processing.  Peak-to-peak 
variation  of  the  original  image  can  be  determined  from  the  degraded 
image  interactively. 

In  our  discussion  of  the  various  methods  of  estimating  the  NSR 
from  the  degraded  image,  we  have  assumed  that  the  noise  process 
is  wide-sense  spatially  stationary  stochastic  process  with  zero  mean 
and  is  spectrally  white,  i.e.,  Sn(u,v)  =  a  constant.  If  F(u,v)  falls 
off  much  faster  in  the  spatial  frequency  domain  than  S^Cu.v),  correctness 
of  the  white  noise  assumption  increases.  We  will  define  the  NSR  as 

NSR  =  (Average  noise  power )/(Average  signal  power)  (131) 

Where  average  signal  power  refers  to  the  average  power  contributed 
by  the  original  image. 
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Estimation  in  the  Spatial  Domain . 


Since  the  noise  and  the  original  image  are  wide-sense  spatially 
stationary  processes,  their  power  density  spectrum  is  defined  by  the 
Fourier  transform  of  their  autocorrelation  function  (Ref.  9),  that 


S  (u,v)  _  /  /  R  (Ax,Ay){exp  -j2Tr(uAx  +  vAy)}  dAx  dAy  (13 

^  _co  _co 


S^(u,v)  A  /  /  R^(Ax,  Ay){exp  -j2rr(uAx  +  vAy)}  dAx  dAy 


It  follows  that  the  autocorrelation  function  is  the  inverse 
Fourier  transform  of  the  power  density  spectrum. 


Rn(Ax,  Ay)  =  ;  /  Sn(u,v)  exp { j 2 -rr( u Ax  +  vAy)}  du  dv 


R^A**  Ay)  =  f  f  S^(u,v)  exp{j27r(uAx  +  vAy)}  du  dv 


From  Equations  (134)  and  (135),  we  can  write 


E {  |h( x , y ) |  }  =  Rn(0,0)  =  f  f  Sn(u,v)  du  dv 

—  OO—OO 


E{  |f(x,y)|  }  =  RpCO.O)  =  /  /  Sc(u,v)  du  dv 


(137 


This  means  that  the  average  noise  power  and  original  image  power 
can  be  identified  by  their  mean  square  value.  Therefore,  we 
may  write  Eq.  (131)  as 


NSR  =  -E(  ln(x>y)l  ^  (138) 

E  {  |f  (x ,  y )  |  / 

since  the  mean  value  of  the  noise  is  zero,  the  numerator  of  the 
Eq.  (138)  is  equal  to 

E{|n(x,y)|2}  =  J^2  (139) 

2 

where  ,  variance  of  the  noise,  can  be  estimated  from  a  relatively 
noisy  region  in  the  degraded  image.  If  a  prototype  image  which  is 
similar  to  that  of  original  image  is  available,  then  we  may  estimate 
the  denominator  term  in  Eq.  (138)  from  the  prototype  image  p(x,y). 

E{  |  f  (x ,  y )  |  2}  =  E{|p(x,y)|2}  (140) 

However,  our  primary  aim  is  to  estimate  all  necessary  parameters  of 
restoration  filter  from  the  degraded  image  itself.  The  degraded 
image  may  be  expressed  as 

g(x,y)  =  fb(x,y)  +  n(x,y)  (141) 

where 

fb(x,y)  =  hd(x,y)  *  f(x,y)  (142) 

If  we  assume  that  the  average  powers  contained  by  fb(x,y)  and  f(x,y)  are 
approximately  equal,  then  we  may  say  that 


E  {  |  f  b(  x ,  y )  |  Z}*  E{|f(x,y)jZ}  (1A3) 

We  can  now  estimate  the  NSR  from  the  degraded  image  itself  as  follows 

E{|f(x,y)|2}  =  E{|g(x,y)|2}  -  E{|n(x,y)|2}  (144) 


where  the  noise  and  the  original  image  are  assumed  to  be  statistically 
independent.  By  using  equations  (141),  (139),  (143),  and  (144) 
we  may  obtain, 


NSR 


E{|g(x,y)|2}-  Cn2 


(145) 


where  the  numerator  and  denominator  denote  the  noise  and  original 

2 

image  average  power,  respectively.  If  O  is  small  compared  to 
2 

E{  |g(x,y)  |  },  then  Eq.  (140)  can  be  approximated  as 

2 

NSR  =  - - - x  (146) 

E{  |g(x ,  y )  |  1 


Estimation  in  the  Frequency  Domain. 

It  is  also  possible  to  estimate  NSR  in  the  spatial  frequency 
domain.  The  power  density  spectrum  of  the  degraded  image  is  given  by 

S  ( u , v )  =  | H^(u, v) | 2  Sf(u,v)  +  Sn(u,v)  (147) 

where  Sn(u,v)  is  constant  for  all  spatial  frequencies  u  and  v  and  is 
approximately  equal  to  the  variance,  i.e.. 
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(148) 


S  (u,v)  %  o  14 

n  n 

By  substituting  Eq.  (148)  into  Eq.  (147)  and  solving  for 
Sn(u, v)/Sf(u, v) ,  we  obtain 


NSR  = 


Sn(u, v) 
Sf(u,v) 


IHd(u,v)l 

- —  a  z 

S  (u,v)-0 
8  n 


(149) 


Since  the  only  constant  in  above  equation  is  O  the  NSR  will  depend 
on  the  spatial  frequencies  u  and  v.  If  the  degrading  OTF  is  zero  at 
any  spatial  frequency,  then  the  NSR  becomes  zero  and  leads  the 
restoration  filter  to  undetermined  ratio.  In  order  to  see  this  effect, 
we  may  insert  Eq.  (149)  into  the  Eq.  (85)  which  the  OTF  of  the 
Wiener  filter. 


Hr(u,v) 


Hd  (u,v) 


lHd(“’V>lV 

+  Sg<u,v)-On2 


(150) 


If  Hd(u,v)  =  0,  then  Hr(u,v)  =  0/0  which  is  undetermined  ratio  that 
we  encountered  in  the  inverse  filtering  method.  We  can  avoid  this  by 
modifying  Eq.  (150)  so  that 


Hr(u,v) 


Hd  (u,v) 


|H  (u,v)|2a2 

|H.(u,v)|2  +-3 - J- 


S  (u,v)-a 
8  n 


if  |Hd(u,v)|  >  R 

if  | H  , (u , v) |  <  R 


where  R=0.  Since  we  cannot  recover  the  zeros  of  the  G(u,v), 

Eq.  (161)  can  be  considered  as  a  reasonable  approximation  to  the 
Wiener  filter.  In  the  absence  of  noise,  it  performs  as  an  inverse 
filter.  Otherwise,  it  deemphasizes  the  noise  component  in  the 
degraded  image. 


IV.  Results 


In  order  to  be  able  to  compare  the  degrading  OTF  estimated  by  one 
of  the  methods  we  discussed  with  the  exact  OTF,  we  may  generate  a  degrading 
OTF  and  convolve  it  with  an  image  to  produce  a  degraded  image  and  then 
estimate  the  degrading  OTF  assuming  a  little  or  no  prior  knowledge 
about  the  degradation. 

To  accomplish  the  above  task,  we  need  to  define  the  spatial  domain 
and  the  spatial  frequency  domain.  We  used  256x256x4-bit  images 
represented  by  a  matrix  of  256  rows  and  256  columns  with  each  element 
of  this  matrix  corresponding  to  a  pixel  in  the  picture  and  each  pixel 
taking  on  gray-scale  values  from  0  (black)  to  15  (white).  In  the 
frequency  domain,  the  discrete  Fourier  transform  of  the  image  is 
represented  by  a  256x256  matrix  of  complex  values.  Both  the  spatial 
(space)  domain  and  the  spatial  frequency  domain  representations  are 
centered  at  the  128'th  row  and  128'th  column  of  the  matrices,  i.e.,  it 
corresponds  to  (x=0,  y=0)  in  the  spatial  domain  and  (u=0,  v=0)  in  the 
spatial  frequency  domain. 

A  simple  computer  graphics  program  was  then  implemented  in 
FORTRAN  5  to  generate  an  image.  This  program  enables  us  to  draw  lines 
from  one  point  to  another  and  circles  centered  at  a  point  in  the  spatial 
domain.  For  example,  the  instructions  below  will  draw  a  line  from 
(x=-50,  y=-50)  to  (x=50,y=50)  with  the  intensity  level  14,  and  a  circle 


at  the  center  (x=0,  y=0)  with  radius  of  10  and  the  intensity  level  7, 
respectively , 

call  LINE( IBUF, 14 , -50,-50, 50,50) 

call  CIRCLE( IBUF ,7,0,0,10) 

where  IBUF  is  an  array  used  to  hold  the  picture  file  (see  Ref.  19). 

A  test  image  was  generated  using  this  program  and  is  shown  in 
Picture  1. 

To  characterize  the  degrading  system  which  we  described  in 
Chapter  II,  "Image  Restoration  System  Model",  we  need  to  generate  a 
point  spread  function  or  an  optical  transfer  function.  We 
implemented  Eq.  (112)  to  obtain  the  0TF  of  the  motion  blur  in  one 
direction.  As  we  can  see  from  Eq.  (112),  it  only  depends  on  the 
frequency  variable  "u"  which  corresponds  to  the  vertical  direction  and 
takes  values  from  -127  to  127  in  the  frequency  domain  that  we  described 
above.  Picture  2  shows  the  0TF  of  the  degrading  system  obtained 
from  Eq.  (112)  for  T=1.0,  and  d=0.021. 

The  degraded  system  is  then  used  to  obtain  a  degraded  image. 

This  degrading  system  simply  convolves  its  point  spread  function 
with  the  input  (original  image)  in  the  spatial  domain.  Since  the 
convolution  in  the  spatial  domain  is  equivalent  to  point-by-point 
multiplication  in  the  spatial  frequency  domain,  the  process  of 
convolving  the  input  (original  image)  picture  with  the  PSF  (point 
spared  function)  to  produce  the  output  (blurred  image)  is  equivalent 
to  first  computing  the  Fourier  transform  of  the  input  (original  image) 
then  multiply  it  with  the  0TF  (optical  transfer  function)  point-by-point. 
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and  then  computing  the  inverse  Fourier  transform  of  the  result  to 
produce  the  output  (.blurred)  image.  This  process  is  summarized  by 
Figure  6  where  2DFFT,  and  2DFFT  1  denote  forward  and  inverse  two- 
dimensional  Fast  Fourier  transform,  respectively  and  "X"  denotes 
point-by-point  multiplication.  Picture  3  shows  the  output  of  this 
degradation  process  for  the  input  image  of  Picture  1. 

re  3 ) 


Figure  5  The  Degrading  (Blurring)  Process 

We  applied  all  filtering  techniques  in  the  Fourier  domain. 

The  input  and  output  picture  files  of  the  two-dimensional  Fast  Fourier 
Transform  software  that  we  used  were  required  to  be  complex  picture 
files,  while  the  input/output  picture  files  of  the  Octek  Image  Analyzer 
system  which  we  used  for  display  of  the  pictures  on  a  monochrome  video 
monitor  were  required  to  be  4-bit  integer  picture  files.  Therefore, 
the  process  of  any  of  the  filtering  operations  that  we  will  discuss 
next  require  the  conversion  of  pictures  from  integer  to  complex  form 
and  vice  versa.  A  general  purpose  filtering  operation  is  summarized 
in  Figure  6  where  VTC  and  CTV  refer  to  the  video-to-cumplex  and 


complex-to-video  conversion  processes,  respectively.  The  filter 
parameters  are  the  NSR  and  the  degrading  OTF. 


Picture  Picture 


Figure  6  Filtering  Operation 

With  our  implementation  on  the  AFIT  Signal  Processing  Lab's  Data 
General  Eclipse  S/250  16-bit  minicomputer:  it  takes  approximately  36 
minutes  (except  for  the  estimation  of  the  parameters  of  the  filter)  to 
process  a  256x256  input  picture  by  a  Wiener  filter  using  above 
algorithm  described  in  Figure  6.  We  will  now  attempt  to  estimate 
the  degrading  OTF  from  the  degraded  image  by  means  of  the  estimation 
methods  described  in  previous  sections. 

As  we  mentioned  in  the  method  of  frequency  domain  inspection, 
we  may  pre-process  the  degraded  image  to  reduce  the  noise  effects. 
However,  we  can  see  by  looking  at  Picture  3  that  there  is  no  apparent 
noise  present.  Therefore,  we  simply  compute  the  Fourier  transform 
G(u,v)  of  the  degraded  image  without  pre-processing.  Picture  4  shows  the 
plot,  of  G(u,v).  The  dark  (black)  points  corresponds  to  the  zeros  of 
the  G(u.v).  Before  examining  the  zeros  of  the  G(u,v),  we  will  turn 
our  attention  back  to  the  degraded  image  (Picture  3).  We  see  that  this 


image  is  blurred  only  in  vertical  direction,  i.e.,  it  contains  a  spread 


in  one  direction  only.  Based  upon  this  observation,  we  will  assume  that 
the  degradation  is  due  to  motion  blur  in  one  direction  which  we  described 
earlier.  If  this  assumption  is  true,  then  we  must  be  able  to  find  a  zero 
pattern  in  the  Fourier  transform  of  Picture  3  which  corresponds  to  the  zeros 
of  the  sine  function  given  by  Eq.  (  90).  Since  the  blur  is  in  vertical 
direction,  we  should  search  for  the  zero  pattern  in  the  vertical 
direction  in  the  plot  of  G(u,v)  shown  by  picture  4,  and  determine  the 
first  zero  of  the  degrading  OTF.  There  are  first  zeros  located  34,  36, 

42,  46,  48,  and  54  pixels  away  (see  Picture  4)  from  the  origin  which  obey 
the  zero  pattern  given  by  Eq .  (90).  However,  only  one  of  these  zeros 
should  be  the  first  zero  of  the  degrading  OTF  given  by  Eq.  (112).  In 
order  to  determine  the  correct  first  zero,  we  may  pick  one  of  those 
zeros,  obtain  the  degrading  OTF  by  Eq.  (112),  and  restore  the  degraded 
image.  If  the  first  zero  that  we  chose  is  not  the  correct  one,  the 
output  of  the  restoration  system  (e.g.,  inverse  filter)  will  not  be  a 
meaningful  picture  and  we  should  repeat  the  process  with  another 
choice  until  we  find  the  zero  that  produces  the  "best"  restoration. 

The  parameter  d  of  the  degrading  OTF  given  by  Eq.  (112)  is  found  by 

d  =  l/|the  first  zero  distance  |  f 152 ) 

A  value  of  d  =  l/4b  gives  the  best  restoration  of  Picture  3  by 
inverse  filtering.  The  degrading  OTF  obtained  by  Eq.  (112)  for  the 
parameter  d  =  l/4b  is  shown  by  Picture  5,  and  the  image  restored 
by  the  inverse  filter  is  shown  in  Picture  6.  When  we  examine  the 


restored  image  (Picture  6),  we  see  noise-like  amplification  artifacts, 
in  addition  to  restored  degradations.  These  artifacts  are  due  to  the 
ill-conditioned  behavior  of  the  inverse  filter.  We  will  consider 
this  problem  later.  When  the  possible  zero  patterns  become  too  many, 
it  takes  too  much  time  and  effort  to  find  the  correct  zero  pattern 
by  this  trial-and-error . 

Since  we  concluded  from  the  degraded  image  (Picture  3)  that  the 
degradation  is  due  to  the  motion  in  vertical  direction,  we  can  employ 
the  line  estimation  method  described  in  Chapter  III,  section  entitled 
"Line  Estimation  of  the  OTF."  If  Equation  (118)  is  modified 
for  motion  in  the  vertical  direction  as 

1  N  i  N 

i  l  CLOG[G(k,  £)  ]  =  CLOG[H  ( H)  ]  +  £  [  CLOG[F(k.  Jl)  ]  (153) 

N  k=l  a  N  k-1 

1  N 

If  the  term  —  £  CLOG[F(k,&)]  is  assumed  to  converge  to  zero,  then 

k-1 

H^il)  is  given  by 

Hd(&)  =  exp  {  i  CLOG[G(k,£)  ]}  k  =  1,  2 . N.  (154) 

As  we  can  see  from  the  Eq.  (154)  Hd(£.)  is  a  column  obtained 
by  averaging  the  columns  of  the  matrix  which  represents  the  complex 
logarithm  of  the  G(u,v),  and  calculating  the  antilogarithm  of  the 
average.  Since  the  degraded  image  is  256x256  in  size,  N  is,  in  this 
case,  equal  to  256  and  H^(£)  is  an  array  (column)  with  256  elements.  To 
estimate  the  first  element  of  the  Hd(£)  corresponding  elements  of  the  256 
columns  of  this  matrix  should  be  summed  and  averaged  and  then 


antilogarithm  of  this  average  should  be  calculated.  This  method 

was  implemented  to  estimate  H^(£)  from  the  complex  logarithm  of  G(u,v). 

Hd(£)  is  then  duplicated  N  times  to  obtain  the  NxN  matrix  of  the 

degrading  OTF.  Picture  7  shows  the  estimated  degrading  OTF.  It  seems 

like  a  sine  function,  but  the  intensity  variations  do  not  obey  a  sine 

1  N 

function  variations  due  to  the  term,  -jr  £  CLOG[F(k,£  )  ] ,  which  we  assumed 

k=l 

to  be  zero.  However,  we  may  obtain  an  exact  sine  function  which  has 
its  first  zero  identical  to  that  of  Picture  7.  This  exact  sine 
function  is  shown  by  picture  8  and  it  is  assumed  to  be  the  degrading  OTF. 
The  degraded  image  is  restored  by  inverse  filter  and  it  is  shown  by 
Picture  9.  We  can  see  the  ill-conditioned  behavior  of  exact  inverse 
filter  in  the  Picture  9.  We  then  modified  the  inverse  filter  as 


follows 


Hr(u,v)  =  A(u, v)  = 


if  Hd(u,v)  >  R 


if  H,(u,v)  <  R 

d 


(155) 


Restored  pictures  using  Eq.  (155)  for  different  R  values  are 


listed  below 


Picture  10  R  =  0.000999 
Picture  11  R  =  0.009999 
Picture  12  R  =  0.099999 
Picture  13  R  =  0.02 


Table  1  Restored  Pictures  and  R  Values  of  Eq .  (155) 


As  we  can  see  from  the  these  restored  pictures,  it  is  possible 
to  overcome  the  ill-conditioned  behavior  of  the  inverse  filter.  For 
small  values  of  R,  the  ill-conditioned  behavior  still  exists,  while  the 
degradation  is  no  longer  restored  for  the  large  values  of  R.  We  could, 
therefore,  find  an  "optimum"  value  of  R  interactively. 

The  degraded  image  (Picture  3)  that  we  restored  in  the  preceding 
experiment  was  a  binary  (two-intensity  level)  picture  generated  by 
computer.  We  next  applied  the  methods  developed  in  Chapter  III  to  the 
degraded  image  shown  in  Picture  14.  This  picture  was  digitized  from 
photographic  film  taken  a  hand  held  35  mm  camera.  The  picture  is  a 
blurred  image  of  the  masking  of  an  aircraft  fuselage.  It  includes  more 
than  two  intensity  levels  with  maximum  and  minimum  intensity  levels, 

14,  0,  respectively.  Its  mean,  meansquare,  and  variance  are 
calculated  to  be  5.11,  34.43,  9.84,  respectively.  Since  no  detailed 
information  was  available  on  how  the  image  was  acquired,  we  assume  no 
prior  information  about  the  degradtion.  Therefore,  we  must  estimate  all 
essential  parameters  from  the  degraded  image  itself.  It  is  not 
possible  to  determine  or  assume  the  form  of  the  degradation  by  looking 
at  the  degraded  image  (Picture  14).  But,  we  can  see  that  it  contains 
noise.  Therefore,  we  will  first  apply  the  Laplacian  operator  to  the 
degraded  image.  Picture  15  shows  the  degraded  image  after  the  Laplacian 
operation.  We  next  take  the  Fourier  transform  of  Picture  15.  The 
magnitude  of  the  Fourier  transform  presented  in  Figure  7  shows  circles  of 
zeros  of  |G(u,v)|  centered  about  the  origin.  The  first  and  second  zero 


are  located  =  3  and  =  6  pixels  away  from  the  origin,  respectively. 
The  ratio  of  R^  to  R2  becomes 

fl  _  3  =  0,  5  (156) 

R9  "  6 

The  degrading  OTF  of  out-of-focus  camera  with  circular  aperture 
forms  approximately  periodic  circular  zero  pattern  and  it  is  given 
by  Eq.  (91).  The  ratio  of  the  first  zero  of  this  pattern  to  its  second 
zero  is  equal  to 


THE  FIRST  ZERO 
THE  SECOND  ZERO 


0.61 

0.117 


0.546106 


(157) 


which  is  approximately  equal  to  the  ratio  R^/R^.  Therefore,  we  will 
assume  that  the  degradation  is  due  to  out-of-focus  camera  with  circular 
aperture  and  its  degraded  OTF  is  of  the  form  of  Eq.  (46).  The  parameter 
r  of  Eq.  (46)  is  obtained  from  the  relationship 

=  3  (158) 

r 

and  then 

r  =  =  0.2033  (159) 

Implementation  of  Eq.  (46)  requires  us  to  calculate  the  Bessel 
function  of  first  kind  of  order  n.  These  series  expansion  of  Jn(x) 
is  given  by  (Ref.  20) 
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Figure  7  The  Zeros  of  |G(u,v) 
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k=0  k!  r( n+k+1 ) 


where  x  =  2irrf,  and  f  =/  u“  +  v“  .  The  spatial  frequency  domain 


that  we  defined  at  the  beginning  of  this  chapter  is  normalized 


by  127  so  that  the  spatial  frequency  variable  can  take  values 


from  -1.0  to  1.0.  This  enables  the  Bessel  function  to  converge  to 


zero.  Eq .  (46)  is  then  calculated  for  N+50,  n=l,  and  r  =  0.2033  to 


obtain  the  degrading  OTF  shown  by  Picture  16. 


After  estimating  the  degrading  OTF,  we  next  need  to  estimate  the 


NSR  in  the  degraded  image.  If  the  NSR  is  defined  by  Eq.  (146),  i.e., 


MSR  = 


E{ I g(x, y )  |  } 


2. 

where  E{|g(x,y)|  }  is  the  mean  square  value  of  the  degraded  image 


(Picture  14)  and  is  computed  to  be  equal  to  34.43,  and  the  variance 


of  the  noise  is  calculated  from  a  relatively  noisy  region  shown  in 


Picture  17  to  be  5.91.  Then,  the  NSR  becomes  0.172.  Picture  18 


shows  the  image  restored  by  the  Wiener  filter  with  degrading  OTF 


of  Picture  16  and  the  NSR  of  0.172.  Another  way  to  obtain  the  degrading 


OTF  is  to  obtain  the  degrading  point  spread  function.  The  Fourier 


transform  of  the  degrading  point  spread  function  gives  us  the 


degrading  OTF.  The  degrading  point  spread  function  given  by 


Eq.  (  45  )  is  calculated  for  R  =  0.203  (shown  in  Picture  19),  then  Fourier 


transformed  to  obtain  the  degrading  OTF  shown  in  Picture  20. 

Picture  21  shows  the  restored  image  obtained  by  Wiener  filter  with  this 
degrading  OTF  (Picture  20)  and  the  same  NSR  (0.172). 

We  also  tested  another  approach  to  estimating  the  degrading  OTF 
presented  in  Eq.  (126).  Since  the  degraded  image  is  256x256,  the 
parameter  N  of  Eq.  (126)  is  equal  to  256.  With  the  degrading  OTF 
estimated  by  Eq.  (126)  and  the  NSR  of  0.172,  the  degraded  image  is  restored 
by  Wiener  filter.  Picture  22  shows  the  restored  image.  As  we  can  see, 
it  is  possible  to  restore  some  degradation,  but  most  of  the  degradation 
is  unrestored.  We  estimated  the  NSR  by  Eq.  (149)  and  used  the  same  OTF 
to  obtain  restoration  filter  OTF  given  by  Wiener  filter.  Restored  images 
with  different  NSR's  are  listed  below 

Picture  23  S  (u,v)  =  5.91 

n 

Picture  24  Sn(u,v)  =  2.4 

Picture  25  Sn(u,v)  =  1.55 

Table  2  Restored  Pictures  and  Noise  Variances 

Where  Sn(u,v)  is  the  noise  power  required  by  Eq .  (149)  to  estimate 
the  NSR.  We  assumed  that  the  noise  is  Gaussian  random  process  with 
zero  mean  and  the  variance  obtained  from  the  noisy  region  (see  Picture  17) 
equal  to  5.91.  The  autocorrelation  of  the  noise  R^  in  the  degraded 
image  is  estimated  by 


noise  power  is  obtained  by 


Sn(u,v)  =|  F{ Rn> 


(162) 


The  NSR  is  then  estimated  by  Eq.  (149).  The  restored  image  is  shown  in 
Picture  26.  In  order  to  avoid  unwanted  artifacts  in  the  restored  images 
with  the  NSR  obtained  from  EQ.  (149),  we  used  modified  Wiener  filter  OTF 
given  by  Eq.  (150)  and  obtained  restored  images  using  different  values 


of  R  values  with  the  same  NSR  value  (see  Pictures  27,  28).  The  NSR 


was  estimated  by  Eq.  (149)  for  S  (u,v)  =  q  =  5.91. 

n  n 


Another  approach  to  estimation  of  the  degrading  OTF  is  the  line 
estimation  method  described  in  the  previous  chapter.  Although  it  was 
developed  for  linear  motion  blur  in  one  direction,  we  will  use  it  for 
estimating  the  degrading  OTF  of  the  degraded  picture.  This  degrading  OTF 
consists  of  only  a  line  either  in  vertical  direction  or  horizontal 
direction.  A  line  in  vertical  direction  can  be  estimated  from  the 


degraded  image  itself  by  Eq.  (120).  In  this  case,  we  have 

,  256 


Hd(k)  =  exp  {  2^  l  CL0G[G(k,£ )  ]  } 

k=l 


(163) 


where  k=l,  2,  3,...,  256.  The  degrading  OTF  obtained  by  Eq .  (120) 
is  shown  in  Picture  29.  The  degraded  image  is  restored  by  Wiener 
filter  with  the  NSR  of  0.172,  and  shown  in  Picture  30.  Secondly,  A  line 
in  horizontal  direction  is  estimated  by 


H.00  =  exp{ 


1 


256 


25  ft 


l  CL0G[G(k ,£ ) ]  } 


(164) 


and  is  shown  in  Picture  31.  The  degraded  image  is  restored  by  Wiener 
filter  with  the  same  NSR(0.172)  and  it  is  shown  in  Picture  32.  As  we 


can  see,  it  is  possible  to  restore  most  of  the  degradation  in  the 

degraded  image.  Therefore,  we  generated  a  sine  function  (Eq.  (112)  with  a 

first  zero  ( d=l / 32 ,  T=1.0)  is  taken  to  the  equal  to  tnat  of  Picture  31. 

Pictures  33  and  34  show  the  sine  function  and  the  restored  image  by 

Wiener  filter  with  the  same  NSR(0.172),  respectively,  Some 

degradation  still  is  present  while  most  of  the  high  frequencies  are 

emphasized  in  the  restored  image  (Picture  34). 

We  finally  tried  to  estimate  the  degrading  OTF  by  segmentation  method 

which  we  described  in  Chatper  III.  The  degraded  image  is  divided  into 

169  subimages  which  are  all  64x64  in  size.  This  is  done  by  forming  a 

window  which  is  64x64  in  size,  then  its  left-hand  corner  placed  at  the  left 

I  hand-corner  (first  row  and  first  column)  of  the  degraded  image,  and 

subimages  are  obtained  shifting  this  window  by  a  shifting  factor  16 

pixels  in  this  particular  case  toward  right  and  down  on  the  degraded 

image.  Each  subimage  is  Fourier  transformed  to  form  the  Eq .  (  93  ). 

Eq.  (96)  is  calculated  by  muliplying  the  magnitude  of  the  Fourier  transfer 

of  each  subimage  point-by-point,  then  taking  1/M  power  of  each  resultant 

^  1/M 

point  and  dividing  by  the  term  [ tt  F.(u,v)]  '  which  is  assumed  to  be 

i 

constant  and  equal  to  1.0.  We  were  not  able  to  restore  the  degraded  imaue 
by  using  the  resulting  OTF.  The  same  algorithm  was  implemented  for  di'Jeren 
constant  values,  but,  it  was  still  not  possible  to  obtain  a  meaningful 
degrading  OTF. 


We  next  estimated  the  power  spectra  given  by  Eq.  (97  )  by  averaging  the 
magnitude  square  of  the  Fourier  transform  of  each  subimage  over  169  sub¬ 
images.  Then  the  power  cepstrum  of  this  estimate  is  obtained.  The 
spikes  occur  periodically  as  in  motion  blur  and  they  are  along  the 
horizontal  direction  (the  degraded  image  (Picture  14)  seems  to  be  blurred 
in  both  horizontal  and  vertical  directions  through).  The  first  spike  is 
4  pixels  away  from  the  origin.  We  assumed  the  degradation  was  due  to 
motion  in  horizontal  direction  and  calculated  the  degrading  OTF  given  by 
Eq .  (ill)  for  d=0.25  it  was  not  possible  to  obtain  a  meaningful 
restoration . 

Another  attempt  to  estimate  the  degrading  OTF  was  made  by  using 
Eq .  (103).  The  degraded  image  was  divided  into  49  subimages  by  shifting 
the  64x64  window  for  32  pixels  in  horizontal  ard  vertical  directions. 

Each  subimage  was  Fourier  transformed  and  the  complex  logarithm  of  the 
Fourier  transform  of  each  subimage  taken  to  form  Eq .  (  99  ) .  The 
iegrading  OTF  was  then  estimated  by  taking  the  average  over  49  subimages 
ind  taking  ant i - logar 1 thm  to  calculate  Eq .  (103)-  The  resultant  degrading 
OTF  seemed  like  a  sine  function  with  its  minimum  at  18  pixel  away  from  the 
center .  We  again  generated  a  sine  function  by  Eq .  (112)  and  used  this 
is  the  degrading  OTF  (shown  in  Picture  35).  and  with  the  same  NSR(0. 172) . 

"he  restored  .rnage  is  shown  in  Picture  36. 

1  be  ibie  to  test  the  methods  that  we  described  in  Chapter  III  for 
*  he  estimation  >t  the  NSR.  we  used  a  V  id  icon  camera  and  the  blurred  image 
•mown  hv  i’:i  tare  The  v  id  icon  camera  was  oriented  to  a  flat  surface  and 


a  picture  was  acquired.  The  light  distribution  on  the  flat  surface 
was  made  as  uniform  as  possible.  Our  aim  was  to  obtain  a  noise  source 
which  might  reflect  the  characteristics  of  the  noise  process  that  we 
discussed  in  Chapter  II  section  entitled  "Noise  Process".  The  noise 
picture  was  then  processed  (Laplacian  operator  applied)  to  emphasize  its 
high  frequency  components  and  subtracted  the  mean.  The  processed  noise 
image  is  shown  in  Picture  37.  Then  the  blurred  image  (Picture  3)  and 
the  noise  image  was  added  together  to  obtain  the  noisy,  degraded  image 
shown  in  Picture  38.  Then  we  used  the  estimation  methods  that  we  describe 
in  the  spatial  domain  and  the  spatial  frequency  domain.  A  noisy  region 
is  selected  as  in  labeled  in  picture.  The  variance  over  this  region  was 
calculated  as  14.7376  and  the  mean  square  calculated  over  the  picture 
as  65.8043.  Then  the  degraded  image  (picture  38)  was  restored  by 
Wiener  filter  using  the  degraded  OTF  shown  in  Picture  5,  and  the 
NSR  estimated  by  Eq.  (103)  as  0.2.  The  restored  picture  is  shown 
in  Picture  40.  As  we  can  see  most  of  the  noise  in  the  degraded  image 
( Picture  38)  is  reduced  in  the  restored  image  (Picture  40)  while  some 
)t  blur  is  left  unrestored.  This  is  because  we  assumed  that  the  power 
-;)**(  t  rum  density  is  white,  i.e.,  constant  and  equal  to  the  variance 
•  •  he  noise  to  derive  a  method  to  estimate  the  NSR  in  Chapter  III. 


is  not  true  in  reality  in  this  case  and  the  restored 


: ma 1 .  Picture  41  shows  the  image  restored  by  power 
.■ar  ion  filter  using  the  same  degrading  OTF  and  the 
i  i  )n-  reduced  noise  in  the  restored  image.  There 


-fir  in  r he  restoration. 


Two  main  approaches  to  estimating  the  leg rad i ng  •  >TF  nave  ?>een 


presented  in  this  thesis.  First,  we  determined  the  lorm  at  •In- 
degrading  OTF  and  whether  it  is  due  to  linear  camera  motion  or  out-ot- 
focus  camera,  then  determined  the  necessarv  parameters  to  estimate  the 
degrading  OTF.  Secondly,  the  degrading  OTF  was  estimated  direct lv  tram  the 
degraded  image.  Methods  to  estimate  the  N’SR  have  been  let  tried  m  spatial 
domain  and  in  spatial  frequency  domain. 

A  trivial  solution  to  the  first  approach  was  the  method  ot 
frequency  domain  inspection.  However,  this  method  required  that  the 
degrading  OTF  must  be  oscillatory  and  the  image  must  contain  appreciable 
high  frequency  components,  well  distributed,  so  that  the  zeros  in  the 
degrading  OTF  can  be  identified.  Furthermore,  the  .VSR  must  be  small 
enough  not  to  abscure  the  zeros;  otherwise,  it  may  become  impossible  to 
determine  its  parameter  (i.e.,  to  recognize  the  zero  pattern)  unless  a 
prior  knowledge  about  the  original  image  is  available.  The  degraded  image 
shown  in  Picture  3  satisfied  these  requirements.  There  were  six  possibility 
for  the  first  zero  of  the  degrading  OTF  that  we  were  searching  for.  we 
had  to  try  each  one  of  them  and  find  out  the  one  which  restores  the  iegradod 
image  best.  The  line  estimation  method  was  then  implemented  to  est  im.it e 
the  degrading  OTF.  We  could  not  estimate  the  exact  degriding  OTF  due 
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wh i .  !i  implement  s  line  estimation  method  while  it  takes  more  than  an  hour 
Mr  the  t  requenc v  domain  inspection  to  ana  1  vze  and  determine  the 
parameters  >t  'he  degrading  OTF.  The  degraded  image  (Picture  3)  was 
•hen  restored  bv  an  inverse  filter.  We  modified  the  inverse  filter  to  avoid 
'he  ill-'  mditioned  behavior.  It  was  possible  to  restore  most  of  the  infor¬ 
mation  which  was  lost  due  to  the  blur.  If  we  compare  the  degraded  image 
i’li  •  ire  l1  »i:h  ‘he  original  image  •  Picture  !/.  most  of  the  information 
;  n  •  .he  iri  les  ind  the  upper  parallel  lines  ot  the  squares  are  lost, 
one  an  no  longer  sav  that  t.hev  are  parallel  lines.  In  the  restored 
pictures,  these  lines  irui  circles,  including  tinv  details,  are 
recovered.  1’his  includes  the  highest  possible  trequencv  data,  the  upper 
parallel  lines  ot  'he  inside  square  which  ire  one  pixel  away  from  each 
alter.  However,  the  restored  images  do  not  look  as  pleasant  and  clear  as 
'he  original  image  does.  Fhev  include  noise  like  variations.  We  had  to 
ompromise  it  this  point,  s i  in  e  we  • ou Id  not  implement  the  exact  inverse 
til' or.  in  'he  ,t  tier  hand,  it  mav  tie  onsidored  more  important 
•  r*"  ivor  ’ho  :  n  t  >rma  f.  :  on  whu  h  :  s  lost  •  nan  •  i  >ht  im  i  more  pleasant. 
im  '  it".  '  irn  o  *  m*  loot  ;  ut  •. -rm.it  ;  on  .  .  r  ivcfi,  ;  •  •  hen  in  he 
or -  -od  »  i  '  '  ':\i'  n ■,  r  :  i :  :i  ’  .  n urn. i: ,  «",•••*  •  :mpt  -ve  '  : . o  ;  u.  i  I  ;  ’  v 


hand  side  of  the  Eq.  (160)  is  theoretically  defined  for  infinite 
number  of  terms  which  we  calculated  it  for  50  terms.  But  this  did  not 
contribute  a  significant  error  to  the  restored  pictures.  The  same 
degrading  OTF  is  calculated  by  the  aid  of  the  Fourier  transform.  The 
degrading  point  spread  function  was  calculated  and  then  Fourier  trans¬ 
formed  to  obtain  the  degrading  OTF.  Some  improvement  was  obtained  in 
the  restoration.  The  segmentation  methods  which  are  successfully  used  to 
restore  the  degraded  images  by  related  references  were  applied  to 
estimate  the  degrading  OTF  from  the  degraded  image.  A  limited 
improvement  in  the  restorations  was  obtained  due  to  the  nature  of  the 
degraded  image.  It  was  possible  to  restore  some  of  the  degradation 
by  the  logarithmic  estimate  method.  This  method  performed  poorly, 
since  it  was  a  rough  approximation  the  degrading  OTF. 

The  estimation  of  the  NSR  in  spatial  domain  performed  well  while 
the  estimation  in  frequency  domain  introduced  an  ill-conditioned 
behavior.  The  ill-conditioned  behavior  of  the  frequency  domain 
estimation  can  be  avoided  by  modifying  the  restoration  filter  OTF 
like  we  did.  Except  the  ill-conditioned  behavior,  the  frequency  domain 
estimation  of  the  NSR  performed  well. 
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V.  Conclusions  and  Recommendations 


Conclusions 

The  method  of  spatial  frequency  domain  inspection  produces  good 
results,  if  the  degraded  image  has  high  frequency  components,  well 
distributed  in  the  spatial  frequency  domain  and  a  relatively  small  NSR 
so  that  the  zeros  of  the  degrading  transfer  function  can  be  identified. 
The  degraded  image  in  Picture  3  met  these  requirements  and  the 
frequency  domain  method  was  successful.  However,  one  must  analyze  the 
zeros  and  consider  many  possible  patterns  for  the  zero  pattern  of  the 
degrading  OTF.  This  can  be  done  for  the  pictures  64x64  in  size,  if  the 
degraded  picture  is  256x256,  512x512,  or  1024x1024,  recognition  of  the 
possible  zero  patterns  may  become  a  very  difficult  and  time  consuming 


task. 

Image  segmentation  methods  required  that  we  segment  the  degraded 
image  into  169  subimages.  It  takes  about  3  hours  of  run-time  for  the 
segmentation  method  to  estimate  the  degrading  OTF  using  the  Signal 
Processing  Lab's  Eclipse  S/250  computer  and  the  program  SEGEST.FR 
given  in  Appendix  B.  We  were  not  able  to  obtain  a  good  estimate  of  the 
degrading  OTF  bv  the  segmentation  methods. 

Line  estimation  method  t  ikes  less  run-time  approx : mat e i v  !| 
minutes)  an«J  et  tort  to  *•••;»  lmat  e  the  legr  iding  '  i'll  .  !  t  v*a>.  lev«>  i  •  .ped 

tor  t  he  teg  r  a  da  t  i  on  due  t  ■  >  in  i  :  ■ .  r  m  line  r  i  m.  a  ;  •  .n  -i  n*  :  -i  .  n  • 


may  be  applied  to  other  pictures  with  unknown  degradations.  If  an 
improvement  is  obtained  in  the  restored  picture,  then  further  steps 
can  be  taken,  such  as  computing  an  exact  sine  function,  or  Bessel  function 
with  the  zeros  of  the  estimated  line.  This  method  performed  well  in  the 
restoration  of  the  degraded  images  that  we  used  in  this  thesis. 

The  logarithmic  estimate  method  performed  poorly,  since  it  was  a 
rough  approximate  to  the  degrading  OTF. 

The  estimation  method  of  the  NSR  in  the  spatial  domain  performed  well 
while  the  estimation  in  the  spatial  frequency  introduced  an  ill-conditioned 
behavior . 

As  a  final  statement  of  this  discussion  we  can  say  that  there  is  a 
chance  to  restore  images  degraded  by  linear  camera  motion  or  out-of¬ 
focus  lens  system  with  circular  aperture  with  no  information  about  the 
degradation  or  the  noise.  The  success  in  restoring  the  degraded  images 
depends  upon  a  good  estimate  of  the  iegrading  OTF  and  the  MSR.  The  success 
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Recommendations 

The  depth  of  this  research  was  necessarily  limited.  Only  limited 
analysis  and  testing  could  be  accomplished  due  to  time  constraints. 
Therefore,  recommendations  for  follow-on  study  consists  of  further 
analysis  and  testing  of  the  estimation  methods. 

A  further  study  may  be  made  on  estimating  the  magnitude  and  the  phase 
of  the  original  image  more  careful Iv  to  obtain  a  better  estimate  ot 
the  degrading  OTF. 

In  the  stages  of  estimating  the  degrading  OTF  bv  segmentation 
method  or  frequency  domain  inspection  method,  the  form  ot  the  degradation 
and  necessary  parameters  to  estimate  the  degrading  OTF  has  '  <>  In- 
obtained  interactively.  A  studv  can  be  made  t  o  reduce  the  human 
i  nterac.  1 1  on  . 
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(Ref.  23,  24).  Therefore,  investigation  and  the  implementation  of 
them  may  be  considered  as  a  further  study. 
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Appendix  B 


Software  Listings  Page 

INVFLT.FR .  B-2 

WIENER.FR .  B-4 

PEQFLT.FR .  B-6 

CONST.FR .  B-8 

SEGEST.FR . B-10 

SEGEST2.FR . B-13 

SEGEST3.FR . B-16 

MTRF.FR . B-19 

DPSF.FR . B-21 

BESSEL.FR . B-23 

LEST.FR  . B-25 

LOGTRF.FR . B-2  7 

CEPST.FR . B-29 

LAP.FR . B-31 

ADV.FR . B-33 

ADC.FR . B-35 

STAT.FR . B-37 

IMAGE.  FR . B-40 

CIRCLE.  FR . B-42 

LINE.  FR . B-43 

SUBIMG.FR . B-45 

SVTC.FR . B-46 

SFFT2.FR . B-47 

UNP.FR . B-47 


B-l 


o  n 


C  < 
C  INVFLT  FR  -  DG  FORTRAN  f  -  LT  MURAT  CINCIOGLU.  JULY  19B5  r. 
C  * 
C  This  program  implements  the  inverse  filter  Inverse  < 
C  Filter  is  modified  to  avoid  the  .ill-conditioned  behaviour  « 
C  due  to  the  leTos  of  the  degrading  OTF  «■ 

< 

The  input  files  will  be  2DFFT  of  the  degraded  image  and  k- 
C  the  degrading  OTF  The  minumum,  value  of  the  degrading  DTF  * 
C  to  avoid  the  leros  will  also  be  asked  The  output  file  «• 
C  will  be  2DFFT  of  the  restored  image  *■ 
C  * 
C  RELOAD  LINE  * 
C  RLDR  INVFLT  ©FLIB0  * 
C  * 
C  COMMAND  LINE  * 
C  INVFLT  * 
C  * 
C  * 


COMPLEX  INI < 1024 ) , IN2< 1024 ) . OUT ( 1024 ) 

INTEGER  INFLNM1 (7) . INFLNM2<7) , 0FLNM1 (7) 

REAL  R 1 j  R2 

C  ACCEPT  INPUT 

ACCEPT "ENTER  BLURRED  IMAGE  Fit  E  NAME  - L " 

READ ( 11,1) INFLNM1 ( 1 ) 

ACCEPT "ENTER  TRANSFER  FUNCTION  FILE  NAME - >" 

READ ( 11,1) INFLNM2 ( 1 ) 

1  FORMAT  <  S 1 3 ) 

ACCEPT "ENTER  OUTPUT  FILENAME  -O  " 

READ ( 11,1 ) OFLNM 1 ( 1 > 

TYPE"  " 

TYPE  "ENTER  MINUMUM  VALUE  TO  AVOID  ZEROS  OF  THE  OTF  " 
ACCEPT"  .  ",R2 

CALL  OPEN ( 1 , INFLNM1. 1, IER) 

CALL  CHECK (IER) 

CALL  OPEN  (2,  INFLNM2,  1,  1LR) 

CALL  CHECK  < IER  > 

CALL  DFILW  ( OFLNM 1 , IER) 

CALL  CFILW  (OFLNMl, 2, KER ) 

CALL  CHECK (KER) 

OPEN  3,  OFLNMl, ATT="OR", ERR-JOO 

100  CONTINUE 
5  TYPE"  " 

TYPE"ENTER  SIZE  OF  FILES  " 

ACCEPT"  256,  120.  64,  OR  32  -  .'".ISJZE 


n 


one 


TEST  FOR  ONLY  T HE  ALLOWABLE  SIZES 

IF ( ISI ZE  EQ  206 )  GO  TO  6 
IF< ISIZE.  EQ  128 j  GO  TO  6 
IF  <  ISIZE.  EQ.  64  )  GO  TO  6 
IF  \ ISIZE  EQ  32)  GO  TO  6 
GO  TO  5 

6  CONTINUE  ;  size  is  okay 

I FACT0R=256/ I S I ZE  ,  compute  scaling  factor 

I T I ME 5=64/ I FAC TOR* *2— 1 

A=0 

DO  3  1=0, ITIMES 

CALL  RDBLK  < 1 ,  16*1,  INI,  16,  IER) 

CALL  CHECK (IER) 

CALL  RDBLK  <  2,  16*1,  IN2,  16,  IER) 

CALL  CHECK  < IER ) 

DO  2  J=1 ,  1024 

R1=REAL( IN2( J) )**2  +  A I MAG  < I N2  <  J ) ) **2 
R 1 =R 1 **0  5 

I F  (  R 1 .  LE.  R2  )  OUT  <  J  )  =  INI  (  J  ) 

IF  <R1  GT  R2  )  OUT (  J  )  =  I N 1 (J)/I N2 ( J ) 

2  CONTINUE 

CALL  WRBLK<3,  16*1,  OUT,  16,  IER) 

CALL  CHECK ( IER  > 

3  CONTINUE 
CALL  RESET 

STOP  "  27T  <7>''7Z  I NVFL  T  FR  " 

END 
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C  ****•*******■***■*•* *•  *  *■*■**■  *•**■*■*•***•■«■  Ir***#^^****####*#***#**#***###**#**#****#*#* * 

^  WIENER.  FR  -  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU,  AUG  1986  * 

C  * 

C  This  program  implements-  the  Wiener  filter  in  the  * 

C  Fourier  domain  The  inputs  are  the  degrading  OTF  and  the  * 

C  Fourier  transform  of  the  degraded  image,  and  the  NSR.  The  * 

C  output  will  be  the  Fourier  transform  of  the  estimated  * 

C  (restored)  image  * 

C  RELOAD  LINE:  * 

C  * 

C  RLDR  WIENER  ©FLIBC  * 

C  * 

C  COMMAND  LINE:  * 

C  * 

C  WIENER  * 

C  Program  will  ask  for  the  degrading  OTF.  degraded  image.  * 

C  and  the  output  file  names  and  the  NSR.  and  the  size  of  the  files.  * 


COMPLEX  INI ( 1024).  IN2< 1024), OUT < 1024) 
COMPLEX  HI, H2, WIENER,  CSN 

t-  INTEGER  INFLNM1 (7).  INFLNM2<7),  0FLNM1 (7) 


ACCEPT  INPUT 

ACCEPT  "ENTER  BLURRED  IMAGE  FILE  NAME- - " 

READ ( 11,1) INFLNM1 ( 1 ) 

ACCEPT" ENTER  TRANSFER  FUNCTION  FILE  NAME - > 

READ( 11,  1 ) INFLNM2  < 1  ) 

1  FORMAT (SI  3) 

ACCEPT"ENTER  OUTPUT  FILENAME 
READ  < 11,1 ) OFLNM 1 ( 1 ) 

TYPE”  " 

ACCEPT"ENTER  N0I5E  TO  SIGNAL  RATIO  =■  ",SN 
CALL  0PEN<1, INFLNM1, 1, IER) 

CALL  CHECK  < IER ) 

CALL  OPEN (2, INFLNM2. 1, ILR) 

CALL  CHECK ( IFR ) 

CALL  DFILW  ( 0FLNM1 ,  IER ) 

CALL  CFILW  (0FLNM1,  2,  KFK) 

CALL  CHECK ( KER ) 

OPEN  3,  OFLNM 1 , ATT  = " OR "  ,  ERP  =  J  00 

100  CONTINUE 
5  TYPE"  " 

,  TYPE"ENTER  SIZE  OF  FILES  " 

ACCEPT"  256,  128,  64,  DR  32 


ISJ  ZE 


i r.,«,:y;Tr.T'">''-'"^yw1 « v  r.i .  r  Fr.'7.7.’7.7  7  V  '■«  '!»  Ml  ^  ^  L'*.1  i'»«ii  u>it< 


TEST  FOR  ONLY  1  HE  ALLOWABLE  SIZES 

IF< ISIZE.  EG  256)  GO  TO  6 
IF ( ISIZE.  EG  120)  GO  TO  6 
IF< ISIZE  EG.  64 )  GO  TO  6 
IF( ISIZE.  EG  32)  GO  TO  6 
GO  TO  5 

CONTINUE  ,  size  is  okay 

I FAC TOR =2 56/ ISIZE 
ITIMES=64/ I FACTOR** 2-1 

DO  3  1=0,  ITIMES 

CALL  RDBLK<  1,  16*1,  INI,  16,  IER  > 
CALL  CHECK  < I ER ) 

CALL  RDBLK  <  2,  16*1.  1N2,  16,  IL'R  ) 
CALL  CHECK (IER) 


compute  scaling  factor 


DO  2  J=l,  1024 

H1=C0NJG< IN2( J) ) 

H2=  IN2< J)*H1 
CSN=CMPLX  <SN,  0  0) 

WIENER=  HI / (  H2  +  CSN  )  , Wiener  Filter  OTF 

OUT (  J  )  =■  I N 1 (J) *W I ENER  ,  2dfft  of  restored  image 

2  CONTINUE 

CALL  WRBLK  <  3,  16*1,  OUT,  16,  IER) 

CALL  CHECK  < I ER ) 

3  CONTINUE 
CALL  RESET 

STOP  "  <:7X7X7>W  I  ENER  " 

END 
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*  c  * 

C  PEQFLT  FR  -  DG  FORTRAN  5  PROGRAM  LT  MURAT  C1NCIOGLU  JUNE  1986  * 
C  * 
C  This  program  implements  the  power  spectrum  equalization  * 
C  filter  Input  files  should  be  the  degrading  OTF  and  the  «- 
C  Fourier  transform  of  the  degraded  image  Output  file  * 
C  will  be  the  estimated  (restored)  image  * 
C  * 
C  RELOAD  LINE:  * 
C  RLDR  PEQFLT  ©FLIBG  * 
C  * 
C  COMMAND  LINE:  * 
C  PEQFLT  * 
C  Program  will  ask  for  input  and  output  file  name,  and  the  NSR.  * 
C  * 
C  * 


C  *********  *r  *****■«•*•*•■*  *•  *■*•#•*-  *■  * 

COMPLEX  INI ( 1024 ),  IN2( 1024).  OUT ( 1024 ) 

COMPLEX  HI. H2. PEQFLT.  CSN 

INTEGER  INFLNM1 (7).  INFLNM2<7>, 0FLNM1 (7) 

ACCEPT  INPUT 

t  ACCEPT-ENTER  DEGRADED  IMAGE  FILE  NAME  — >" 

READ ( 11,1) INFLNM1 ( 1 ) 

ACCEPT "ENTER  OTF  FILE  NAME  — 7 " 

READ ( 11,1) I NFLNM2  < 1  ) 

1  FORMAT (SI  3) 

ACCEPT "ENTER  OUTPUT  FILE  NAML  — >" 

READ ( 11,1 )0FLNM1 ( 1 ) 

TYPE"  " 

ACCEPT-ENTER  NOISE  TO  SIGNAL  RATIO  =  ",SN 
CALL  OPEN ( 1 , INFLNM1 ,1,1 FR ) 

CALL  CHECK ( IER) 

CALL  OPEN <2, INFLNN2, 1- IER) 

CALL  CHECK (IER) 

CALL  DFILW  (0FLNM1 ,  IER) 

CALL  CFILW  (0FLNM1 , 2, KOR ) 

CALL  CHECK < KER ) 

OPEN  3,  0FLNM1,  ATT="OR",  ERR-100 

100  CONTINUE 
5  TYPE"  " 

TYPE"ENTER  SIZE  OF  FILES  ” 

ACCEPT"  256,  128,  64,  OR  32  — >”.  ISIZE 


B-6 


non 


TE5T  FOR  ONLY  THE  ALLOWABLE  SIZES 

IF<  ISIZE.  EQ.  256)  GO  TO  6 
I F  (  ISIZE.  EQ.  128)  GO  TO  6 
IF<  ISIZE.  EQ.  64)  GO  TO  6 
IF  (  ISIZE.  EQ.  32)  GO  TO  6 
GO  TO  5 

6  CONTINUE  i  size  is  okay 

IFACT0R=256/ ISI ZE  ;  compute  scaling 

ITIMES=64/IFACT0R**2-1 

DO  3  1=0,  IT I ME S 

CALL  RDBLK ( 1 ,  16*1,  INI,  16,  I ER ) 

CALL  CHECK < IER ) 

CALL  RDBLK  ( 2,  16*1,  IN2,  16,  IER) 

CALL  CHECK < IER ) 

DO  2  J-l,  1024 

HI =CONJG ( IN2  <  J )  ) 

H2=  I N2  <  J ) *H 1 
CSN=CMPLX (SN,  0  0) 

PEQFLT=  <  < 1  0/ (H2+CSN) ) ) **0  5 

OUT ( J ) *PE«FLT 

2  CONTINUE 

CALL  WRBLK<3,  16*1,  OUT,  16,  IER ) 

CALL  CHECK  < IER ) 

3  CONTINUE 


factor 


CALL  RESET 

STOP " <7><7><7>PEGFLT  FR 


o  r>  o 


1 '  •  ' »  W  Wm1 v j^jv  A  WWW  v  w 


TTjrmwrqrwCTTi 


C a-*******************************-#*****-***************'*** **#**•*•*•#■*<  te* 


C  v 
C  CONST  FR  -  FORTRAN  5  PROGRAM  LT  MURAT  C  INC  I  OGLU  MAY  B6  t- 
C  *- 
C  The  NSR  is  estimated  in  the  spatial  frequency  domain.  * 
C  Wiener  filter  is  then  used  to  restore  the  degraded  image.  < 
C  Inputs  are  2DFFT  of  the  degraded  iinsge.  and  the  noise  power  * 
C  Output  u/ill  be  the  restored  image  * 
C  * 
C  RELOAD  LINE:  * 
C  RLDR  CONST  ©FLIB©  * 
C  * 
C  COMMAND  LINE:  * 
C  CONST  * 
C  Program  will  ask  for  input  and  output  file  name/  create  the  * 
C  output  file  and  then  ask  for  file  size.  * 
C  * 


C*****  ************•■*••**■**•#■•*•****■**■#••#■*•#•*•  *******************************.•**■ 

COMPLEX  G(  1024),  F<  1024),  HW,  HBPS.  NSR,  DENUM,  GPS,  NPS,  SIZE 
I NTEGER  I NFLNM 1 <  7 ) ,  OFLNM 1 <  7 ) ,  I NFLNM2 (  7 ) 

COMPLEX  HB ( 1024 ) 

REAL  NLEV.R.M 

ACCEPT  INPUT 

ACCEPT" ENTER  DEGRADED  IMAGE  FILE  NAME  — >" 

READ (11, 1 ) I NFLNM 1 < 1 ) 

ACCEPT "ENTER  TRANSFER  FUNCTION  FILE  NAME  — >" 

READ <11,1)  INFLNM2  < 1 ) 

1  FORMAT ( SI  3 ) 

ACCEPT"ENTER  OUTPUT  FILENAME  ->” 

READ (11, 1 ) OFLNM 1 < 1 ) 

CALL  0PEN(1, I NFLNM 1 , 1, IER) 

CALL  CHECK (IER) 

CALL  OPEN <2, INFLNM2, 1. IER) 

CALL  CHECK (IER) 

CALL  DFILW  ( 0FLNM1 , IER) 

CALL  CFILW  (0FLNM1, 2, KER) 

CALL  CHECK (KER) 

OPEN  3.  0FLNM1,  ATT="OR", ERR=100 

100  CONTINUE 
5  TYPE"  " 

TYPE"ENTER  SIZE  OF  FILES  " 

ACCEPT"  256,  128,  64,  OR  32  — ISIZE 

TYPE"  " 

ACCEPT"  ENTER  MINUMUM  OF  DEG  OTF  >",R 

ACCEPT"  ENTER  NOISE  POWER  ... NLEV 


TEST  FOR  ONLY  THE  ALLOWABLE  SIZES 


IF<  ISIZE.  EQ.  256)  GO  TO  6 
IF  < ISIZE  EQ  12B )  GO  TO  6 
IF  (  ISIZE.  EQ.  64)  GO  TO  6 
IF<  ISIZE.  EQ.  32)  GO  TO  6 
GO  TO  5 

CONTINUE  ;  size  is  okay 

IFACT0R=256/ ISI ZE  , compute  scaling 

ITIMES=64/IFACT0R**2-1 

SI ZE=CMPLX  <  ISIZE**2,  0  0) 

NPS  =  CMPLX  (NLEV,  0.  0) 

DO  3  1=0, ITIMES 

CALL  RDBLK<  1,  16*1,  G,  16,  IER  ) 

CALL  CHECK (IER) 

CALL  RDBLK  <  2,  16*1, HB,  16,  IER) 

CALL  CHECK (IER) 

DO  2  J=l, 1024 

CALCULATE  NOISE  TO  SIGNAL  RATIO  ( NSR ) 

HBPS  =  HB(J)*C0NJG<  HB ( J )  )  *NPS 

GPS  =G(U)*CONJG<  G(J)  ) 

DENUM  =  GPS  -NPS 
NSR  =HBPS/DENUM 

CALCULTE  RESTORATION  FILTER  TRF 
M=  HB ( J  >  *CONUG ( HB  <  J ) ) 

M=  SQRT(M) 

HW  =  CONJG ( HB ( J ) ) /  (  CONJG ( HB ( J ) ) *HB ( J )  +  NSR 
IF  (  M.  LT.  R  )  HW  =  CMPLX  <1  0,0  0 ) 

CALCULATE  RESTORED  IMAGE 

F(J>=  G ( J )  *  HW 

2  CONTINUE 

CALL  WRBLK  <3,1 6* I ,  F,  16,  IER  ) 

CALL  CHECK (IER) 

3  CONTINUE 
CALL  RESET 

STOP  "  <7><;7><:7>C0NST.  FR  “ 

END 


factor 


o  o 


C  ***** 

C 

C 


***•**»*******•«•**■«■*  ******■»****«■**«■  ****■»#****#**##****#********•«■***«•***♦# 

« 

SEGEST.FR  -  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  19B6 


,y 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


This  program  implements  the  segmentation  method 
Its  input  is  the  degraded  image  It  divides  the  image  into 
square  regions  (subimages)#  each  subimage  is  then  converted  to 
complex  representation#  and  Fourier  transformed  Finally  thr 
magnitude  squares  of  the  Fourier  transform  of  each  sub  image  are 
summed  and  averaged  by  the  number  of  subimages 
LOAD  LINE 

RLDR  SEGEST  SUB  IMG  SVTC  SFFT2  FFT2  F5PICBUF  LB  @FLIB@ 


COMMAND  LINE: 


* 

* 

* 

* 

* 


SEGEST 

Program  will  ask  for  the  I/O  file  names  and  two  temporary 
file  names#  size  of  the  subimages,  and  the  shifting  factor  of  the 
window  which  subimages  are  obtained  by  shifting  this  window 
The  output  file  name  will  be  the  estimated  power  spectra 


* 

* 

* 

* 

* 

* 


SUMS- 


**»*#********#*#*«*****#**#**t*****#***#****t*******«*************** 


INTEGER  DIMGNM  <  7 ) .  TMPNM1 (7),  TMPNM2(7>. OTFNM ( 7 ) 
INTEGER  R. C. DELTHA,  M.  ISTR,  ISTP, OPTION 
COMPLEX  INI (1024).  IN2(1024),  OUT (1024) 

COMPLEX  CM, CMAG 

> 

9  I/O  CONSTANTS 

PARAMETER  KTTYIN=1 1 
PARAMETER  KTTYOUT^IO 
C 

ACCEPT "ENTER  DEGRADED  IMAGE  FILE  NAME 
READ(KTTYIN, 40)  DIMGNM ( 1 ) 

ACCEPT "ENTER  OUTPU i  FILE  NAME: .  >" 

READ(KTTYIN, 40)  OTFNM ( 1 ) 

ACCEPT "ENTER  TEMPORARY  FILE  NAME  #1  ....>" 
READ(KTTYIN. 40)  TMPNM 1(1) 

ACCEPT "ENTER  TEMPORARY  FILE  NAME  #2.  .  .  .  >" 

READ ( KTTY IN, 40)TMPNM2( 1 > 

40  FORMAT (S40) 


10  TYPE "ENTER  SIZE  OF  THE  SUB  IMAGES" 


ACCEPT"  128.64,32  .  .  >".  ISIZE 

ACCEPT "ENTER  SHIFTING  FACTOR  DELTHA 

IF ( ISI ZE. EQ  128)  GO  TO  20 

IF( ISIZE  EQ  64)  GO  TO  20 

IF( ISIZE  EQ.  32)  GO  TO  20 

GO  TO  10 

CONTINUE 


20 


random  file 


CALL  DFILUKOTFNM, IER) 

CALL  CFILW<  OIPNM,  2, IER) 

CALL  CHECK (IER) 

OPEN  7,  OTFNM,  ATT="OR ",  ERR- 120 
CONTINUE 

CALL  CF ILW < TMPNM 1 . 2, IER)  ;  initialize 
CALL  CHECK (IER) 

CALL  CF  ILW  (  TMPNM2,  2,  IER) 

CALL  CHECK (IER) 

ISTR-1  ,  starting  pointtl, 1)  for  first  segment 

I STP=256- I SIZE 
DELTHA=DELTHA- 1 

M=0  ; number  of  segments 

OPTION=l  ,  2D  forward  Fourier  transform 

I FAC TOR =2 56/ ISI ZE 
ITIMES=64/IFACT0R**2-1 

DO  2  1=0- IT IMES 
DO  3  J=l. 1024 

OUT ( J ) =CMPLX ( 0  O- 0  O) 

CONTINUE 

CALL  WRBLK ( 7 .  16*1-  OUT-  16-  IER) 

CALL  CHECK (IER) 

CONTINUE 

DO  25  R=l-  ISTP,  DELTHA 

DO  35  C  =  l,  ISTP,  DELTHA 
M=M+1 


CALL  SUBIMG(DIMGNM,  TMPNM1,  ISIZE, C,  R)  i  segmentation 

CALL  SVTC< ISIZE. TMPNM1, TMPNM2)  ;  convert  to  complex 
CALL  SFFT2(TMPNM2- OPTION)  ;  2D~foruiard  FFT 

OPEN  2,  TMPNM2,  ATT= "OR " , ERR= 1 0 1 
OPEN  7,  OTFNM,  ATT= "ORR " , ERR= 1 01 
CONTINUE 

DO  55  1=0,  I TIMES 

CALL  RDBLK  <  2,  16*1,  INI,  16,  IER) 

CALL  CHECK (IER) 

CALL  RDBLK  ( 7,  16*1,  IN2,  16,  IER) 

CALL  CHECK (IER) 


DO  50  J=l, 1024 


OUT  <  U)=IN2<  J)  +  INI  <  J  >  *CON JG  (  INI  <  J )  ) 
CONTINUE 

CALL  WRBLK(7,  16*1,  OUT,  IER) 

CALL  CHECK (IER) 

CONTINUE 

TYPE"SEGMENT  #",  M 
CONTINUE 


CONTINUE 

CM=CMPLX (FLOAT (M),  0  0) 

DO  77  1=0,  ITIMES 

CALL  RDBLK(7, 16*1, OUT, IER ) 

CALL  CHECK (IER) 

DO  87  J=l,  1024 


OUT ( J ) =OUT ( J ) /CM 


;  average  by  li 


CONTINUE 

CALL  WRBLK(7, 16*1, OUT, IER) 
CALL  CHECK (IER) 

CONTINUE 

TYPE" NUMBER  OF  SEGMENTS12",  M 
CALL  RESET 

STOP " <7><7><7>SEGEST  FR " 

END 


v  * .  •.  \  \  - 


c ************************************************************** ************** 

c  * 

V.  SEGEST2.  FR  —  DG  FORTRAN  5  -  LT  MURAT  Cl  NCI  OGLU.  AUG  19Q6  * 

£ 

C  This  program  implements  the  segmentation  method.  * 

C  Its  input  is  the  degraded  image  It  divides  the  degraded  image  * 

C  into  M  subimages.  The  magnitude  of  the  Fourier  transform  of  each  * 

C  subimage  is  then  multipled  and  1/M  power  of  the  result  is  * 

C  calculated  to  estimate  the  magnitude  of  the  degrading  OTF.  * 

C  RELOAD  LINE:  * 

C  * 

C  RLDR  SEGEST2  SVTC  SFFT2  FFT2  UNP  F5PICBUF.  LB  ©FLIB©  * 

C  * 

C  COMMAND  LINE:  * 

C  * 

C  SEGEST2  * 

C  Program  will  ask  for  the  I/O  and  two  temporary  * 

C  file  names,  the  size  of  the  subimages  and  the  shifting  factor.  * 

C  The  output  will  be  the  magnitude  of  the  estimated  degrading  OIF.  * 


C *********************************** ***************************************** 


m 


c 


40 

C 

C 

10 


20 


INTEGER  DIMGNM  <  7 ) ,  TMPNM1 (7), TMPNM2<7>. OTFNM ( 7 ) 
INTEGER  R. C, DELTHA. M. ISTR, ISTP, OPTION 
COMPLEX  INI < 1024),  IN2< 1024), OUT( 1024) 

REAL  Rl.  R2,  RM 
I/O  CONSTANTS 
PARAMETER  KTTYIN=1J 
PARAMETER  KTTYOUT^IO 

ACCEPT" ENTER  DEGRADED  IMAGE  FILE  NAME  :....>" 
READ<KTTYIN. 40)  DIMGNM( 1 ) 

ACCEPT "ENTER  OUTPUT  FILE  NAME. .  >" 

READ<KTTYIN. 40)  OTFNM  < 1 ) 

ACCEPT"ENTER  TEMPORARY  FILE  NAME  41  ....>" 
READ(KTTYIN, 40)  TMPNM1 < 1 ) 

ACCEPT "ENTER  TEMPORARY  FILE  NAME  #2.  .  .  >" 

READ ( KTTYIN,  40)TMPNN2< 1 ) 

FORMAT (S40) 


TYPE"ENTER  SIZE  OF  THE  SUB  IMAGES" 

ACCEPT"  128,64,32  .  >", ISIZE 

ACCEPT "ENTER  SHIFTING  FACTOR  :  ...  .  DELTHA 

IF<  ISIZE.  EQ.  128)  GO  TO  20 

IF< ISIZE  EQ.  64)  GO  TO  20 

IF<  ISIZE.  EQ.  32)  GO  TO  20 

GO  TO  10 

CONTINUE 


B—  1 3 


CALL  DFILW(OTFNM, I EH > 

CALL  CF ILW  (  OTFNM,  2,  IER)  ;  2  random  file 

CALL  CHECK (IER) 

OPEN  7,  OTFNM,  ATT="OR”, ERR=120 
CONTINUE 

CALL  CFILW( TMPNM1 , 2,  IER)  i  initialize 
CALL  CHECK < IER) 

CALL  CFILW ( TMPNM2,  2,  IER) 

CALL  CHECK ( IER ) 

ISTR=1  i  starting  point(l,l)  for  first  segment 

ISTP=256-ISIZE 

DELTHA=DELTHA- 1 

M=0  inumber  of  segments 

OPTION=l  ;  2D  forward  Fourier  transform 

I FACTOR =256/ ISI ZE 
ITIMES=64/IFACT0R**2-1 

DO  2  1=0, I TIMES 
DO  3  J=l, 1024 

OUT  <  J ) =CMPL  X  < 1 .  0, 0.  0)  i  initialize 

CONTINUE 

CALL  WRBLK(7.  16*1, OUT,  16,  IER) 

CALL  CHECK (IER) 

CONTINUE 

DO  25  R=l. ISTP, DELTHA 

DO  35  C=l,  ISTP,  DELTHA 
M=M+1 

CALL  SUB  I MG  <  DIMGNM,  TMPNM1 ,  I  SI ZE, C, R  )  .segmentation 

CALL  SVTC ( ISIZE, TMPNM1,  TMPNM2)  ;  convert  to  complex 
CALL  SFFT2CTMPNM2.  OPTION)  ;  2D-forward  FFT 

OPEN  2, TMPNM2,  ATT=" OR ” , ERR= 1 0 1 
OPEN  7,  OTFNM,  ATT= "ORR ” ,  ERR= 1 01 
CONTINUE 

DO  55  1=0. ITIMES 

CALL  RDBLK  ( 2,  16*1,  JN1,  16,  IER) 

CALL  CHECK (IER) 

CALL  RDBLK ( 7 ,  16*1.  IN2-  16,  IER) 

CALL  CHECK (IER) 


DO  50  J=l,  1024 

Rl  =  <  REAL  < I N 1 ( J ) )  )**2  +  (  AIMAG(INl(J 
R2= (  REAL ( I N2 ( J ) )  >**2  +  <  AIMAG(IN2(J 
R  1  =  R1**0  5 
R2=  R2**0.  5 
RR=  R 1 *R2 

OUT ( J ) =CMPLX ( RR ,  0  0) 

CONTINUE 

CALL  WRBLK<7, 16*1, OUT, IER ) 

CALL  CHECK (IER) 

CONTINUE 

TYPE  "SEGMENT  #'\M 
CONTINUE 


CONTINUE 

RM= 1 .  0/ (FLOAT <M) ) 

DO  77  1=0, ITIMES 

CALL  RDBLM7,  16*1,  OUT,  IER) 
CALL  CHECK (IER) 

DO  87  J=l, 1024 

RR=(  REAL < OUT (J))  )**  RM 

□UT( J)=CMPLX (RR,  0.  O) 

CONTINUE 

CALL  WRBLK(7, 16*1, OUT, IER) 
CALL  CHECK (IER) 

CONTINUE 

TYPE-NUMBER  OF  SEGMENT5=",  M 
CALL  RESET 

STOP " <7><7><7>SEGEST2. FR " 


)  )  ** 

)  )  ** 


nj  nj 


n  o 


C  *********************************************  #■*■*■*■*••*•***•*•***-*#■*•***•***■#••#■•#••*>■■#•■»*•**■*  f 


SEGEST3.  FR  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  19B6 


s 


This  program  obtains  implements  the  segmentation  * 
method  Obtains  an  estimate  of  the  complex  logarithm  of  degrading  * 
OTF  by  avereging  the  complex  logarithm  of  the  Fourier  transform  of  * 
subimages  * 


RELOAD  LINE: 


RLDR  SEGEST3  SUB  IMG  SVTC  UNP  SFFT2  F5PICBUF.  LB  ®FLIB@ 


COMMAND  LINE: 


SEGEST3  * 

Program  will  ask  for  I/O  file  and  two  temporary  file  names  * 
and  for  the  shifting  foctor.  and  the  size  of  the  subimages  * 


C  ******■»*«■■«■*  ************************  ****************************  *****'M-* *--*■»<••***■* 


INTEGER  DIMGNM ( 7 ) . TMPNM1 (7),  TMPNM2<7).  OTFNM ( 7 ) 
INTEGER  R,  C. DELTHA, M. ISTR, ISTP, OPTION 
COMPLEX  INI < 1024), IN2< 1024), OUT( 1024) 

COMPLEX  CM 

REAL  ROUT,  IOUT, RR,  II 


I/O  CONSTANTS 
PARAMETER  KTTYIN=1 1 
PARAMETER  KTTY0UT=10 


ACCEPT" ENTER  DEGRADED  IMAGE  FILE  NAME 
READ(KTTYIN, 40)  DIMGNM ( 1 ) 

ACCEPT" ENTER  OUTPUT  FILE  NAME: . >" 

READ(KTTYIN,  40)  OTFNM  < 1 ) 

ACCEPT "ENTER  TEMPORARY  FILE  NAME  #1 
READ<KTTYIN, 40)  TMPNM1 ( 1 ) 

ACCEPT "ENTER  TEMPORARY  FILE  NAME  #2.  .  .  >" 

READ ( KTTYIN, 40)TMPNM2< 1 ) 

FORMAT (S40) 


vlv  20 


TYPE "ENTER  SIZE  OF  THE  SUB  I MAGES' 
ACCEPT"  128,64,32 

ACCEPT"ENTER  SHIFTING  FACTOR 
IF<  ISIZE.  EG  128)  GO  TO  20 
IF<  ISIZE.  EG.  64)  GO  TO  20 
IF<  ISIZE  EG.  32)  GO  TO  20 
GO  TO  10 
CONTINUE 


. >", ISIZE 
>",  DELTHA 


V  'j  ■>  v  7  7’;  vw/  w/,7.?.’Vivjy.  y.  *.'  JL1  r,1  e.1  * ;  *  ■. » ■„.»  n.»  \ 7  «.■».■  v  v 


c 

CALL  DFILW ( OTFNM,  IER) 

CALL  CFILW  ( OTFNM,  2, IER)  i  2  random  file 
CALL  CHECK (IER) 

OPEN  7,  OTFNM,  ATT="OR",  ERR=120 
120  CONTINUE 


CALL  CF ILW < TMPNM 1 , 2, IER)  ;  initial 
CALL  CHECK (IER) 

CALL  CF ILW ( TMPNM2, 2, IER) 

CALL  CHECK (IER) 

ISTR  =  1  i  starting  point(l,  1)  for  first  segment 

ISTP=256- ISIZE 
DELTHA=DELTHA- 1 

M=0  inumber  of  segments 

0PTI0N=1  i  2D  forward  Fourier  transform 

I FAC TOR =256/ I  SI ZE 
ITIMES=64/IFACT0R**2-1 


DO  2  1=0.  ITIMES 
DO  3  J=l, 1024 

OUT  <  J )  =CMPLX  (0.0,00) 

3  CONTINUE 


CALL  WRBLK  (  7,  16*1,  OUT,  16,  IER) 
CALL  CHECK (IER) 

CONTINUE 

DO  25  R=l, ISTP, DELTHA 

DO  35  C=l. ISTP, DELTHA 
M=M+1 


99  CALL  SUB  I MG ( DIMGNM, TMPNM1 ,  ISIZE, C, R)  i  segmentation 

CALL  SVTC( ISIZE. TMPNM1, TMPNM2)  ;  convert  to  complex 
CALL  SFFT2(TMPNM2,  OPTION)  ,  2D-foru»ard  FFT 

OPEN  2, TMPNM2,  ATT="QR”, ERR=101 
OPEN  7,  OTFNM,  ATT="ORR”, ERR=101 
101  CONTINUE 


DO  55  1=0, ITIMES 

CALL  RDBLK ( 2, 16*1, INI, 16, IER) 

CALL  CHECK (IER) 


CALL  RDBLK  (  7,  16*1,  IN2,  16,  IER) 
CALL  CHECK (IER) 


r*i 


DO  50  J=l,  1 024 


OUT  <  J )  =  IN2 ( J )  +  CLOG  ( I N1 ( J ) ) 

50  CONTINUE 

CALL  WRBLK(7, 16*1. OUT, IER) 
CALL  CHECK (IER) 

55  CONTINUE 

TYPE" SEGMENT  #",M 
35  CONTINUE 


25  CONTINUE 

109  CM=CMPLX  ( FLOAT  <  M  ),  0.  0) 

DO  77  1=0,  ITIMES 

CALL  RDBLK ( 7,  16*1 , OUT ,  IER ) 

CALL  CHECK (IER) 

DO  87  J=l,  1024 


OUT( J)=OUT( J) /CM  ; average  by  M 

CONTINUE 

CALL  WRBLK(7, 16*1, OUT, IER) 

CALL  CHECK (IER) 

CONTINUE 

TYPE-NUMBER  OF  SEGMENTS= ",  M 
CALL  RESET 

STOP " <7><7><7>SEGES I  3.  FR" 

END 


o  o  o  o 


*********#*#***»*#*#-***************************************************#*#** 

* 

MTRF  FR  -  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  1986  * 

* 

This  program  obtains  a  degrading  OTF  for  the  motion  * 

rC  blur  in  the  horizontal  direction.  The  degrading  OTF  of  the  * 

C  motion  blur  in  vertical  direction  can  be  obtained  by  tranforming  * 


C  the  output  file.  * 
C  * 
C  RELOAD  LINE:  * 
C  * 
C  RLDR  MTRF  @FLIB©  * 
C  * 
C  COMMAND  LINE.  * 
C  * 
C  MTRF  * 
C  Program  will  ask  for  the  degrading  OTF  file  name.  * 
C  total  exposure  time  <T).  and  total  displacement  distance  (d).  * 
C  and  the  size  of  the  output  file.  * 
C  * 
C  * 


C  ******************************  ******  ****  ****•*#*##■***•*•**  #***■**•***■*•*■#■•#■#■****##* 


COMPLEX  TRF  <  256 ) 

PARAMETER  KTTYIN=1 1. KTTY0UT=10,  PI  =  3.  1415927 
REAL  Rl.  R2,  T,  d.  R,  RO.  X 
INTEGER  OFLNM ( 7 ) 


WRITE  <  KTTYOUT ,  12) 

12  FORMAT  (  "ENTER  FILE  NAME . >'\Z) 

READ<KTTYIN.  14)  OFLNM < 1 ) 

14  FORMAT <  S39  ) 

TYPE"  " 

ACCEPT "ENTER  EXPOSURE  TIME . >".  T 

TYPE"  " 


ACCEPT  "ENTER  TOTAL  DIPLACEMENT  <"d").  .  .  >'\  d 

CALL  DFILWCOFLNM. IER ) 

CALL  CFILW<OFLNM, 2, HER) 

CALL  CHECK (  KER  ) 

CALL  F0PEN<2. OFLNM) 

C  CALCULATE  TRANSFER  FUNCTION 

DO  20  U*l.  256 
X=  U-128  0 
R0=  PI*X*d 
IF (  X  EQ.  0.  0  )  R=T 

IF  (  X  NE  0.  0  )  R=<  T  *S  I N  (  RO  )  )/(  RO  ) 

R 1  =  R*COS (  RO  ) 

R2=  -R*SIN<  RO  ) 

TRF(U)  =  CMPLX  <  Rl. R2  ) 

20  CONTINUE 


TYPEWRITING  TRANSFER  FUNCTION" 
DO  30  1=1/ 256 

DO  40  J=l> 256 
WR ITE  <  2 )  TRF(J) 
CONTINUE 

IF  (  MOD<  I,  50).  EG.  0)  TYPE  "DONE"/ 
CONTINUE 
CLOSE  2 

SI  OP " <C7><7><C7>MTRP .  FR“ 

END 


oo>o  uit-*  non 


C***#*#HHt**##**»#*******#***#*«*#*******#*#******************###*<-** 


C  * 
C  DPSF  FR-  DG  FORTRAN  5  -  LT  MURAT  Cl  NCI  OGLU-  MAY  1986  * 
C  * 
C  This  program  calculates  two-d  i  men  s  i  oana  1  point  spread  function* 
C  for  an  out-of-focus  camera  with  circular  aperture  in  the  * 
C  normalized  coordinates.  * 
C  * 
C  * 
C  RELOAD  LINE:  * 
C  RLDR  DPSF  @FLIB©  * 
C  * 
C  COMMAND  LINE:  * 
C  DPSF  * 
C  Program  will  ask  for  an  output  file  name,  create  the  * 
C  output  file  and  then  ask  for  file  size,  the  radius  of  the  * 
C  defocusing  error.  * 


C  ******  **■«•**■•#■  *******#*************■*****•*•**•»•**•#••#■■*•«••#•*•«•*•*•*■#••»(•*•*•»••*•*  fr  *  * 


COMPLEX  OUT (1024) 
INTEGER  0FLNM1 < 7 ) 

REAL  R.  RD.  VALi  VAL2 
REAL  X,  Y,  RSI  ZE,  OR 
PARAMETER  PI=3.  1415927 


ACCEPT  INPUT 
1  FORMAT (SI  3) 

ACCEPT "ENTER  OUTPUT  FILENAME 
READ (11, 1 ) 0FLNM1 < 1 ) 

CALL  DFILW  ( 0FLNM1 , IER ) 

CALL  CFILW  ( 0FLNM1 , 2, KER ) 

CALL  CHECK (KER) 

OPEN  3,  0FLNM1, ATT="OR", ERR=100 

00  CONTINUE 
TYPE"  " 

TYPE  "ENTER  SIZE  OF  FIIES  " 

ACCEPT"  256,  128,  64,  OR  32  — ISIZE 

ACCEPT "RAD I US  OF  DEFOCUSING  — >" > R 


TEST  FOR  ONLY  THE  ALLOWABLE  SIZES 

IF< ISIZE  EQ  256)  GO  TO  6 
IF ( ISI ZE.  EQ  128)  GO  TO  6 
IF ( ISIZE. EQ  64 )  GO  TO  6 
IF( ISIZE  EQ.  32)  GO  TO  6 
GO  TO  5 

6  CONTINUE  ;  size  is  okay 

RSI ZE  =  ISIZE  i  type  conversion 


B— 2 1 


C *********************************** ******************* ************** ******** 

c  * 

BESSEL. FR  -  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  1986  * 

C  This  program  obtains  a  degrading  OTF  for  defocused  * 

C  camera  with  circular  aperture.  * 

C  * 

C  RELOAD  LINE:  * 

C  * 

C  RLDR  BESSEL  ®FLIB€  * 

C  * 

C  COMMAND  LINE:  * 

C  * 

C  BESSEL  * 

C  Program  will  ask  for  the  degrading  OTF  file  name.  * 

C  radii  of  the  first  minima,  number  of  terms  to  approximate  * 

C  the  Bessel  function  of  first  kind  with  order  one.  and  the  * 

C  size  of  the  output  file.  * 

C  * 

C ******************************************************************** ******** 

COMPLEX  TRF ( 256 ) 

PARAMETER  KTTYIN=1 1 . KTTYOUT=10 

REAL  C<256)  ,  X .  R.  F,  U.  V,  OR.  C2  (  256 ) .  WX,  WY.  SI  ZE.  APRX ,  J1 .  J.  N.  K 
INTEGER  OFLNM < 7 ) 

REAL  K2  .  A 
C 

flc  WRITE  (KTTYOUT,  12) 

fb  FORMAT ( "ENTER  FILE  NAME  FOR  THE  TRANSFER  FUNCTION  >",Z> 

READ(KTTYIN.  14)  OFLNM' 1) 

14  FORMAT <  S39  ) 


TY  E"  " 

TV PE" ENTER  DEFOCUSING  ERROR (  Radii  of  the  first  minima)" 
ACCEPT"  >",R 

TYPE"  " 

TYPE "HOW  MANY  TERMS  SHOULD  BE  CALCULATED  FOR  THE  BESSEL  FUNC 
ACCEPT"  > " , APR  X 


ACCEPT "ENTER  SIZE  OF  THE  OUTPUT  FILE 


>".  SIZ 


CALL  DFILW(OFLNM.  IER) 
CALL  CFILWCOFLNM,  2,  HER) 
CALL  CHECK (KER) 

CALL  F0PEN<2.  OFLNM) 


CALCULATE  TRANSFER  FUNCTION 

OR=SI ZE/2  0 

DO  50  1  =  1 , APRX 

C  < I )  =  1 . o 

C2< I )=1  0 

CONTINUE 

DO  100  1=1-  APRX 
X=  I 

DO  200  J=l.  0-  X 
C ( I ) ®C ( I ) * J 
C2< I ) =C2 ( I ) *U 
CONTINUE 

C2( I ) =C2  < I )*<  X  +  l.  0) 

CONTINUE 
N=  APRX 

DO  20  WY=1  O,  SI ZE 
DO  30  WX= 1 .  0, SIZE 

U=  < WY— OR ) /OR  ;  normalized  coord. 

U= ( WX— OR ) /OR  i translate  org.  to  the  cntr.  of  the  file 

F= (  U**2+V**2  )  **0.  5 
X=  2.  0*R*F 

IF  <  X .  EQ.  0.  0)  GO  TO  310 
Jl  =  X/2  0 
A=l.  0 

DO  300  1=1. APRX 

K=  I  ;  type  conversion 

K2=  1 .  0+2.  0*K 

IF<MOD<  I.  2).  EG.  0)  A=-1.0 
Jl~  Jl-A*  <<X/2.  0>*»K2  >  / <C< I >*C2( I > ) 

CONTINUE 

IF(X.  EQ.  O.  0)  TRF(WX  )=CMPLX  (  1.  0-0  0) 

IF  ( X  EQ  0.0)  GO  TO  19 
J1=J1/(R«F) 

IF  (  X.  NE.  0.  O)  TRT  (  WX  )  =  CMP  LX  <  J 1 . 0.  O) 

WRITE (2)  TRF(WX) 

CONTINUE 

IF(MOD(WY,  OR)  EG  O  0)  TYPE "DONE.  .  :"» WY 
CONTINUE 
CLOSE  2 

STOP ” <7><7><7>BESSEL  FR " 

END 


V  V  V  1  *  \«  v  '  w’-  «  y**>  • 


'C  LEST  FK - DG  FORTRAN  5  -  L.T  MURAT  CINCIOGLU.  AUG  1986  + 

C  * 

C  This  prooTam  implements  the  line  estimation  method  * 

C  It  estimates  a  line  in  the  vertical  direction  and  duplicates  * 

C  this  line  to  obtain  an  output  file  which  is  identical  in  size  * 

C  to  the  input  file  * 

C  * 

C  RELOAD  LINE.  * 

C  * 

C  RLDR  LEST  <?FL  I  BE  * 

C  * 

C  COMMAND  LINE:  * 

C  * 

C  LEST  * 

C  Program  will  ask  for  the  input  file  size,  and  input/  * 

C  output  file  names  * 

C  * 

c  * 

c  * 


c  **********<•****■«■**«■**«■**«■*  **■**■«■■«■****■**■*■#•*#*#********•******#**#***  ********  *** 


COMPLEX  INI <256)  .  OUT  1 <256). OUT <256) 

COMPLEX  NB. EST 

INTEGER  INFLNM <  7 ) ,  OFLNM  <  7 ) , UMP.  K,  M,  N,  STEP,  WINDOW.  KSIZE 

C  ACCEPT  INPUT 

TYPE" INPUT  FILE  SHOULD  BE  FOURIER  TRANSFORM  OF  " 

TYPE" THE  DEGRADED  IMAGE  .  " 

TYPE"  " 

WR  I TE  <  10,  12) 

12  FORMAT< "ENTER  INPUT  FILE  NAME  =  — >",Z) 

READ  < 11,  14)  I NF  LNM  < 1  ) 

14  FORMAT <S39) 

WRITE< 10, 16) 

16  FORMAT < "ENTER  OUTPUT  FILE  NAME  =  — >".Z) 

READ  < 11,  14  >  OFLNM < 1  ) 

TYPE"  ENTER  SIZE  OF  INPUT  FILE  : " 

ACCEPT"32,  64,  123,  256  — >",  ISIZE 

TYPE"  " 

C  TYPE"  ENTER  WINDOW  SIZE  " 

C  ACCEPT"1, 2,  4,  8,  16,  32.  64,  128  T  •• ,  WINDOW 

C 

C 


WINDOW 


1 


WRITE (2)  OUT ( J ) 

NTINUE 

(  MOD (1.64)  EQ  0 ) 7  YPE "  DONE 


ODD 


c  *  ******■►*■«■**•#■**  ***•#•**  ♦  ■fr-**-***********-*-*  *  «■ 


c  * 

C  LOGTRF  FR  DG  FORTRAN  5-  LT  MURAT  CINCIOGLU,  JUNE  85  * 
C  * 
C  This  program  estimates  degrading  OTK  from  degraded  *• 
C  image  and  calculates  Wiener  filter  transfer  function  * 
C  for  the  restoration  filter  * 
C  * 
C  LOGTRF  estimates  the  OTF  by  taking  complex  logarithm  of  * 
C  Fourier  transform  of  the  degraded  image  and  deviding  by  * 
C  total  number  of  pixels  in  a  line.  «- 
C  * 
C  RELOAD  LINE:  * 
C  RLDR  LOGTRF  SFLIBG  * 
C  < 
C  COMMAND  LINE:  * 
C  LOGTRF  * 
C  Program  will  ask  for  input  and  output  file  name,  create  the  * 
C  output  file  and  then  ask  for  file  size.  * 
C  Input  file  name  should  be  Fourier  transform  of  the  degraded  * 
C  image.  No i se-to-Si gna 1  ratio  is  also  asked  for  Wiener  filter.  * 
C  No i se-to-s i gna 1  raio  is  assumed  to  be  constant.  * 
C  * 


COMPLEX  HB  < 1024 ) ,  G  < 1 024 ) , HW ( 1 024 ) . OUT ( 1 024 ) 
INTEGER  I NFLNM1 <  7 ) ,  0FLNM1 ( 7 ) 

REAL  NS 

ACCEPT  INPUT 

ACCEPT "ENTER  DEGRADED  IMAGE  FILE  NAME  ....>" 
READ ( 11.1) INFLNM1 ( 1 ) 

1  FORMAT  <  S 1 3 ) 

ACCEPT" ENTER  OUTPUT  FILENAME  ->" 

READ ( 11,1 ) 0FLNM1 ( 1 ) 

CALL  0PEN(1, INFLNM1 <  1, IER) 

CALL  CHECK (IER) 

CALL  DFILW  ( 0FLNM1 , IER ) 

CALL  CFILW  ( 0FLNM1 . 2. KFR ) 

CALL  CHECK (KER) 

OPEN  3,  0FLNM1 , ATT="OR " >  ERR = 1 00 

100  CONTINUE 
5  TYPE"  " 

TYPE "ENTER  SIZE  OF  FILES  " 

ACCEPT"  256,  128,  64,  OR  32  — >", 


ISIZE 


TEST  FOR  ONLY  THE  ALLOWABLE  SIZES 


IF<  ISIZE.  EQ.  256)  GO  TO  6 
IF<  ISIZE.  EQ.  12B)  GO  TO  6 
IF(  ISIZE.  EQ.  64)  GO  TO  6 
IF <  ISIZE.  EQ.  32)  GO  TO  6 
GO  TO  5 

CONTINUE  ;  size  is  okay 

IFACT0R=256/ 1 SI ZE  ;  compute  scaling  factor 

ITIMES=64/IFACT0R**2-1 

SIZE  =  CMP LX ( ISI ZE**2,  0.0) 

CNS  =  CMPLX  (NS,  0.  0) 

DO  3  1=0, ITIMES 

CALL  RDBLK(  1,  16*1,  G,  16,  IER  ) 

CALL  CHECK (IER) 

DO  2  J=l, 1024 

HB(J)  =  CLOG(  G(J)  )  /  SIZE  i  estimated  OTF  of  degradi 
HW(J)  =  CONUG (  HB  <  J )  )/(  HB ( J ) *CONJG ( HB ( J ) )  +  CNS  )  ;  W 
OUT ( j )=  HW ( j ) *G  <  J  ) 

2  CONTINUE 

CALL  WRBLK  <  3,  16*1,  OUT,  16,  IER) 

CALL  CHECK (IER) 

3  CONTINUE 
CALL  RESET 

STOP " <7><7><7>L0GTRF " 


non  o'  ^  o  o 


C  * 
C  CEPST. FR  —  DG  FORTRAN  5  PROGRAM  LT  MURAT  CINCIOGLU  5  JUNE  1986  * 
C  i. 
C  This  program  computes  the  cepstrum  of  the  degraded  image  * 

C  The  input  file  should  he  the  Fourier  transform  of  the  * 
C  degraded  image.  *: 
C  * 
C  RELOAD  LINE.  *• 
C  RLDR  CEPST  SFFT2  ©FLIBC  * 
C  * 
C  COMMAND  LINE:  * 
C  CEPST.  FR  * 
C  Program  mill  ask  for  input  and  output  file  name*  create  the  * 
C  output  file  and  then  ask  For  file  size  < 
C  * 


C  **********  *■  ****************•■*•**■»<■■*******-***•****  ■&■******•*****■*•»■*■**  •#•*■«•■«■  r  *  * 

COMPLEX  INI ( 1024), OUT < 1024 ) 

COMPLEX  CEP 
REAL  AMP 

INTEGER  INFLNM1 <  7 ) , OF LNM 1(7) 

ACCEPT  INPUT 

ACCEPT  "ENTER  INPUT  FILENAME  #1 
READ  < 11.1) INFLNM1 ( 1 ) 

1  FORMAT  <  S 1 3 ) 

ACCEPT "ENTER  OUTPUT  FILENAME 
READ  < 11,1 )0FLNM1 ( 1 ) 

CALL  0PEN<1, INFLNM1, 1, IER) 

CALL  CHECK ( IER ) 

CALL  DFILW  ( 0FLNM1 , IER) 

CALL  CFILW  ( 0FLNM1 , 2, KER ) 

CALL  CHECK < KER) 

OPEN  8,  0FLNM1 , ATT= "OR " , ERR=  J  00 

00  CONTINUE 
TYPE"  " 

7 YPE "ENTER  SIZE  OF  FILES  . 

ACCEPT"  256,  128,  64,  OR  32  —  O " ,  ISIZE 

TEST  FOR  ONLY  THE  ALLOWABLE  SIZE5 

IF ( I S I ZE  EQ  256)  GO  TO  6 
IF(  ISIZE.  EQ.  128)  GO  TO  6 
IF  < ISIZE  EQ.  64)  GO  TO  6 
I F  <  ISIZE  EQ.  32)  GO  TO  6 
GO  TO  5 

6  CONTINUE  ;  size  is  okau 


f  /  / 


■  V.V  V  V  V  V" 


J=l>  ISIZE 


DO  4  K=l, ISIZE 
READ  < 1 )  INI <  K ) 

AMP  =  REAL  < INI <  K ) )  +  A I  MAG ( INI <  K ) ) 

AMP=SQRT( AMP ) 

AMP=LOG< AMP ) 

IF<  AMP  NE.  0.  0)  AMF‘-  LOG  ( AMP  ) 

IF  (AMP.  EQ.  0.  0)  AMP  =  0.  0 
CEP=CMPLX ( AMP,  0.  0) 

OUT  <  K ) =CEP 

WRITE(B)  OUT ( K ) 

CONTINUE 

CONTINUE 
CLOSE  8 

CALL  SFFT2 ( 0FLNM1 , 1 ) 


STOP " <7><7><7>CEPST.  FR' 
END 


_>  *  .'jTji  'j  m  '  J  _  . 


*  ws.*' •  v -.■V/ vV.- vV/VvS  A 


"  \‘V-  ■ 


C************************-**  +  ********«-*<r****«**fe*««*«*«**«*0-*  *-*****-«*«  ««**«-** 

c 

C  LAP  FR  —  PG  FORTRAN  0  L T  MURAT  CINCIOGLU.  AUG  1986 

C 

C  This  prograffi  applies^  the  Laplacian  operator  to 

C  the  input  image 

C 

C  RELOAD  LINE:  ^ 

C 

C  RLDR  LAP  F5PICBUF. LB  @FL IB©  • 


COMMAND  LINE: 


C  *****■****■&•***■#■***••*■**•**■*•*•*•  <-  * 


PARAMETER  NBUFSZ=300,  KTTY0UT=10,  KTTY I N= 1 1 
PARAMETER  NC0LSMAX=256 
INTEGER  I BUF ( NBUFSZ ) 

INTEGER  I  ARRAY ( NCOLSMAX ,  3 ) ,  I ARRAYS  (  NCOLSMAX ) ,  INAME<20) 


INITIALIZE  INPUT  AND  OUTPUT  BUFFERS  AND  READ  IN  INPUT  PICTURE 
WRITE  (KTTYOUT, 996) 

FORMAT < "INPUT  PICTURE  FILE  NAME  =  " ,  Z ) 

READ  (KTTYIN, 997)  I NAME ( 1 ) 

FORMAT (S39) 

CALL  PICFMT< IBUF, INAME, IFLC) 

CALL  GFMT ( IBUF,  NROWS, NCOLS,  NB ITS,  1 MODE ) 

IF  ( IMODE  EQ  3)  GO  TO  10 

IF  < (NROWS  LT.  3) .  OR.  (NCOLS  LT  3)  )  GO  TO  10 
IF  (NCOLS  GT  NCOLSMAX)  GO  TO  10 

I F  ( NB ITS.  GT.  11)  GO  TO  10  ,  to  preclude  integer  overflow 
CALL  MAKB(IBUF) 

CALL  P  ICIN( IBUF,  INAME,  IFLG ) 

INITIALIZE  CONVOLUTION  PROCESS 

CALL  GROW  < IBUF,  1,1, NCOLS,  I  ARRAY ( 1,1)) 

CALL  GROW ( IBUF,  2,  1, NCOLS,  I  ARRAY (1,2)) 

11=0 

12=1 

13=2 

NR  1 =NROWS- 1 
NC1=NC0LS-1 
I ARRAY2 ( 1 ) =0 
I ARRAY2  <  NCOL S )=0 


B— 3 1 


PERFORM  CONVOLUTION  A  ROW  AT  A  TIME 


DO  200  1=2,  NR  1 

IF  <MOD<  I,  10)  EQ.  0)  WR I TE ( KTTYOUT,  99B  >  I 
FORMAT ( “ROW  =  “,13) 

I  1 =MQD  < II,  3 )  +  1 
1 2=M0D ( 12,  3) +  1 
1 3=M0D  < 13, 3) +  1 

CALL  GROW  (  1 BUF ,  I  1 ,  1,  NCOLS,  I  ARRAY  (  1,  13)  ) 

DO  lOO  J=2,  NCI 

IVAL=  I  ARRAY  <  J-l ,  I  1  )  *0  -I  ARRAY  (  J,  1 1  )  +  IARRAY  (  J+l ,  1 1  ) 

I VAL=IVAL-IARRAY  <  J-l ,  12 )  + 1 ARRAY ( J,  12 ) *5- I ARRAY < J+i ,  12) 

I VAL=I VAL+I ARRAY < J-l ,  1 3 ) *0- 1  ARRAY ( J,  1 3 )  +  1 ARRAY < J+ 1 ,  I3)*0 
ifdVAL.  LT  0)  IVAL=<-1  )*IVAL 
if(IVAL.  GT.  15)  I VAL=  1  5 

IARRAY 2(J)=  IVAL  ,  (laplacian  of  image) 

CONTINUE 

CALL  PROW < I BUF ,  I,  1, NCOLS,  IARRAY2) 

CONTINUE 

OUTPUT  RESULT  AND  RELEASE  BUFFERS 
WRITE  (KTTYOUT, 999) 

FORMAT ( “OUTPUT  PICTURE  FILE  NAME  =  ”,Z) 

READ  (KTTYIN. 997)  I NAME ( 1 ) 

CALL  P ICOUT < I BUF, I NAME, IFLG ) 

CALL  RELEU  IBUF) 

STOP " <7><7><7>LAP .  FR " 

END 


of  image) 


•W  ■_*  •,'ji  v  vs>  ^ 


c  * 

~  ADV.FR  --  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU,  AUG  1986  * 
C  * 
C  This  program  adds  two  integer  picture  files  to  * 
C  obtain  the  output  file  * 
C  * 
C  RELOAD  LINE:  * 
C  * 
C  RLDR  ADV  F5PICBUF.  LB  ©FLIB©  * 


COMMAND  LINE. 


C  ***************•**■**■*****■■*■**■#■■«•*■*■<!■****■  *  *♦***■»**•»•***********•#•  ■#•**•#■•* •*•*•*•*****•**#* * 


PARAMETER  KTTYIN=1 1, KTTY0UT=10, NBUFSZ=300 
PARAMETER  MAXCOL  =-256,  MAXR0W=256 

INTEGER  IBUF1 (NBUFSZ).  I BUF2 <  NBUFSZ ) *  IBUF(NBUFSZ) 
INTEGER  I NAME 1 (7), INAME2(7>, ONAME ( 7 ) 

INTEGER  I ARRAY 1 (MAXCOL ) , I ARRAY2( MAXCOL ) . I ARRAY ( MAXCOL ) 


WRITE ( KTTYOUT , 990) 

FORMAT ( IX, "ENTER  PICTURE  FILE  NAME  #1 
READ<KTTYIN, 991 )  I NAME 1 < 1 ) 

FORMAT <S39) 


>",  Z) 


WRITE (KTTYOUT,  993) 

FORMAT ( IX, "ENTER  PICTURE  FILENAME  #2 
READ(KTTYIN, 991 )  INAME?(1) 


WRITE (KTTYOUT, 992) 

FORMAT < IX, "ENTER  OUTPUT  FILE  NAME 
READ(KTTYIN, 991 )  ONAME ( 1 > 


— >  ",  1) 


CALL  P I CFMT ( I BUF 1 ,  I NAME  1 ,  I HDR ) 

CALL  MAKB(IBUFl) 

CALL  PICIN( IBUF1 , INAME1 , IHDR ) 

CALL  GFMT ( IBUF1,  NROWS, NCOLS,  NB ITS,  IMODE) 


CALL  P I CFMT < I BUF2- INAME2, IHDR) 

CALL  MAKB ( I BUF2 ) 

CALL  P I C I N ( I BUF2,  I  NAME 2,  I HDR ) 

CALL  GFMT ( I BUF2, NROWS, NCOLS,  NBITS,  IMODE) 
CALL  PFMT ( I BUF ,  NROWS,  NCOLS,  NBITS,  IMODE) 
CALL  MAKB (I BUF) 


DO  500  1=1, NROWS 

CALL  GROW < I DUF 1 ,  1 ,  1 ,  NCOLS,  I ARRAY 1 )  ;  get  rows 

CALL  GROWC IBUF2, I, 1, NCOLS, IARRAY2) 

DO  600  J=  1,256 

IVAL  = I ARRAY 1 ( J )  +  I  ARRAY? ( J )  ;  add  rows 

IF  < I VAL.  LT  0 )  IVAL  =-l*IVAL 
IF(  IVAL.  GT.  15)  I VAL=  1  5 
I ARRAY  <  J) =1 VAL 
CONTINUE 

I F  <  MOD  < 1 , 64 ) .  EQ  0 )  TYPE "DONE  .  :",I  ;  report  to  keyboard 

CALL  PROW( IDUF,  I,  1,  NCOLS,  IARRAY) 

CONTINUE 

CALL  PICOUT < IDUF, ONAME, 0) 

CALL  RELD  < IBUF1 ) 

CALL  RELB  < IBUF2 ) 

CALL  RELB  < IBUF ) 

CALL  RESET 

STOP " <7><7><7>ADV " 


non 


C  *******#*■*■#*#*#  **•*■*  *•*■•*★■**•*•  **■#■*■*■*****#■*■*•**■*■■#■***•*•*«■*.#•*#■»**#•**#•*•****#-.*<.  * 


c  * 

C  ADC.  FR  -  DC  FORTRAN  5  -  LT  MURAT  Cl NCI OGLU  AUG  1986  * 
C  * 
C  This  program  adds  two  complex  image  to  obtain  the  * 
C  output  file.  * 
C  * 
C  RELOAD  LINE:  * 
C  RLDR  VCM  ©FLIBe  * 
C  « 
C  COMMAND  LINE:  * 
C  ADC  * 
C  Program  will  ask  for  input  and  output  file  name,  create  the  * 
C  output  file  and  then  ask  for  file  size.  * 
C  * 


C  **************  ************************  **************-**************-K-  ■* 

COMPLEX  INK  1024),  IN2< 1024) ,  OUT < 1024 ) 

INTEGER  INFLNM1 ( 7 ) ,  INFLNM2 (7 ) ,  0FLNM1 (7) 

ACCEPT  INPUT 

ACCEPT” ENTER  INPUT  FILENAME  #1  ->" 

READ( 11, 1 ) INFLNM1 ( 1 ) 

ACCEPT-ENTER  INPUT  FILENAME  42  ->" 

READ  <11,  1 ) INFLNM2  < 1 ) 

1  FORMAT (SI 3) 

ACCEPT "ENTER  OUTPUT  FILENAME 
READ ( 11,1 )0FLNM1 < 1 ) 

CALL  0PEN(1, INFLNM1 , 1, IER ) 

CALL  CHECK (IER) 

CALL  OPEN (2, INFLNM2, 1, IER) 

CALL  CHECK (IER) 

CALL  DFILW  ( 0FLNM1 , IER ) 

CALL  CFILW  ( OFLNM  1,2,  KER  ) 

CALL  CHECK (KER) 

OPEN  3,  OFLNM 1 , ATT=”OR " , ERR= 1 00 

100  CONTINUE 
5  TYPE"  " 

TYPE"ENTER  SIZE  OF  FILES: " 

ACCEPT"  256,  128,  64,  OR  32  — >",  ISI ZE 


C  TEST  FOR  ONLY  THE  ALLOWABLE  SIZES 
C 

IF<  ISIZE.  EQ.  256)  GO  TO  6 
IF<  ISIZE.  EQ.  128)  GO  TO  6 
IF(  ISIZE.  EQ.  64)  GO  TO  6 
IF<  ISIZE.  EQ.  32)  GO  TO  6 
GO  TO  5 

6  CONTINUE  i  size  is  okay 

IFACT0R=256/ ISI ZE  ; compute  scaling  factor 

ITIMES=64/IFACT0R**2-1 

DO  3  1=0, ITIMES 

CALL  RDBLK  < 1 ,  16*1,  INI,  16.  IER  > 

CALL  CHECK (IER) 

CALL  RDBLK  <  2,  16*1,  1N2,  16,  IER) 

CALL  CHECK (IER) 

DO  2  J=l, 1024 

OUT ( J )  =  I N 1 ( J )  + 1 N2  <  J ) 

2  CONTINUE 


CALL  WRBLK ( 3, 16*1, OUT, 16, IER) 
CALL  CHECK (IER) 

3  CONTINUE 


CALL  RESET 

STOP " <7><7><7>ADC " 

END 
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C ************************************** ************************************** 

c  * 

C  STAT  FR  -  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  1986  * 

f-S  * 

•w  This  program  obtains  the  mean,  var  iance,  mean  squat r,  * 

C  maximum,  and  minumum  values  of  a  picture  or  a  given  region  * 

C  in  the  picture  * 

C  RELOAD  LINE:  * 

C  * 

C  RLDR  STAT  F5PICBUF  LB  ©FLIB©  * 

C  * 

C  COMMAND  LINE:  * 

C  * 

C  STAT  * 

C  Program  will  ask  for  the  picture  file  name  and  * 

C  the  c  oord  inates  <  r  oui,  c  o  1  )  of  the  region  where  the  mean,  * 

C  variance,  mean  square  to  be  calculated.  * 


C ***************************************** *********************************** 


REAL  MEAN,  MEANSQ,  VAR,  STDEV,  PR,  SN,  L(  16) ,  SV,  N 
INTEGER  K,  V,  Y,  M,  MIN,  MAX 

PARAMETER  NBUFSZ=300,  KTTY0UT=10,  KTTYIN=11 
PARAMETER  NC0LSMAX=256 
INTEGER  I BUF ( NBUFSZ ) 

INTEGER  I NAME (20) 


C  INITIALIZE  INPUT  BUFFER  AND  READ  IN  PICTURE  FILE 


10  WR I TE  <  KTT YOUT ,  996 ) 

996  FORMAT < "ENTER  INPUT  PICTURE  FILE  NAME  =",Z) 
READ<KTTYIN, 997)  I NAME ( 1 ) 

997  FORMAT (S39) 

CALL  PICFMT( IBUF, INAME, IFLG) 

CALL  GFMT  < IBUF,  NROWS, NCOLS, NBITS,  IMODE) 

IF<  IMODE.  EQ.  3)  GO  TO  10 

IF( (NROWS  LT.  3)  OR  (NCOLS  LT.  3))  GO  TO  10 
IF < NCOLS  GT.  NCOLSMAX)  GO  TO  10 
I F  ( NBITS  GT.  1  1  )  GO  TO  10 
CALL  MAKB(IBUF) 

CALL  PICIN( IBUF, INAME, IFLG) 


TYPE"  ENTER  REGION  BOUNDARIES" 

ACCEPT "ROW  BOUNDARY  FROM — >".K 
ACCEPT"  TO  — :'".M 
ACCEPT "COLUMN  BOUNDARY  FROM — } "  .  V 
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MAX=0.  0 
MI N=  1  5.  O 
N=0  0 


DO  2  1=0, 15  ;  initialize 

L ( I ) =0  0 

2  CONTINUE 


DO  20  I=K,  M  .  count  possible  intensity  levels 

DO  78  U=V,  Y 

CALL  GPNT ( IBUF,  I»  J,  I VAL ) 

IF ( I UAL. LT  MIN)  MIN=IVAL  ;  max, min  search. 

IF< IVAL.  GT.  MAX )  MAX=IVAL 

IFdVAL.  EG  0)  L<0)=L  <0)4-1.  0 
IF  <  IVAL.  EG  1  )  L<1)=L  <1)4-1.  0 
IF<  IVAL.  EG.  2)  L<2)=L  <2)4-1.  0 
IF  <  IVAL.  EG  3 )  L  <  3 )  =L  <3)4-1.  0 
IF<  IVAL.  EG.  4)  L(4)=L  <4)4-1.  0 
IF  <  IVAL  EG.  5 )  L<5)=L<5)4-1.  0 
IF<  IVAL.  EG  6)  L<6)=L<6)4-1.  0 
IF  <  IVAL.  EG  7 )  L<7)=L  <7)4-1.  O 
[)  IF(IVAL  EG  8)  L<8)=L<8)4-1.  0 

IF<  IVAL.  EG.  9)  L<9)=L<9)4-1  0 
IF  (IVAL  EG  10)  L<  10)=L<  10)4-1  0 
IF  <  IVAL.  EG.  11  )  L<  11  >=L<  11  >4-1.  0 
IF  <  IVAL  EG  12)  L<  12)=L<  12)4-1.  0 
IF  (IVAL.  EG  13)  L<13)=L<  13)4-1.  0 
IF(IVAL  EG  14)  L<  14)=L<  14)4-1.  0 
IFdVAL  EG  15)  L<15)=L<  15)4-1.  0 
N=N4i  o  ;  number  of  elements 

78  CONTINUE 

20  CONTINUE 

MEAN=0  O 
VAR=0  0 
MEANSQ=0.  0 


88 


DO  88  1=0  0,  15  O 
PR  =L<I)/N 
MEAN=  MEAN4-  I*PR 
VAR  =  VAP-K  < I -MEAN ) *#2  ) *PR 
MEANSQ=MEANSQ4-  (  1**2  )*HR 

CONTINUE 


STDEV=VAR**0  5 
SN=VAR /MEANSQ 


REPORT  MEAN,  MEAN  SQUARE,  VARIANCE,  MAX.,  MIN. 
TYPE"  " 


TYPE"MEAN  =",  MEAN 
TYPE"  " 

TYPE”VARI ANCE  =",VAR 


TYPE"  " 

TYPE “STANDART  DEVIATION=", STDEV 
TYPE"  " 

TYPEMMEANSQUARE  =" ,  MEANSG 
TYPE"  " 

TYPE " ( VAR  I ANCE ) / < MEANSQUARE )  =  " , SN 
TYPE"  " 

TYPE"MINIMUM=", MIN 
T YPE ” M A  X I MUM= " , MAX 

STOP " <7><7><7>ST AT . FR" 

END 


o  o  o  o  o 


C*********  ***************************  <-r********r**fr*****-****-fr*****-p-fr********< 

c  * 

C  IMAGE.FR  —  DG  FORTRAN  5  -  LT  MURAT  CINCIOGLU.  AUG  1986  * 

C  * 

C  This  program  is  used  to  create  the  test  image  (Picture  1)  * 

C  * 

C  RELOAD  LINE:  * 

C  * 

C  RLDR  IMAGE  LINE  CIRCLE  F5PICBUF.  LB  @FLIB@  * 


COMMAND  LINE: 

IMAGE 

Program  mill  ask  for  the  output  file  name  and  the 
background  intensity  level. 


C  ********************  ********  *************************************** ********* 


PARAMETER  KTTY I N“ 1 1 » KTTYOUT^IO, N3UFSZ=300 
INTEGER  IBUF ( NBUFSZ ) * IARRAY(256) 

INTEGER  XI,  X2.  Yl, Y2,  I NAME  <  7 ) ,  IBGL,  IVAL 


CALLS  SUBROUTINE  L I NE  <  I  BUF,  I  VAL,  X  1 ,  Y 1 ,  X2,  Y2 ) 
CALLS  SUBROUTINE  C I RCLE ( IBUF ,  I VAL,  X , Y,  RADI OUS ) 

X1CX2.  Y1CY2 


wr ite (KTTYOUT, 12) 

format  <  "ENTER  FILE  NAME . >”,Z) 

read (KTTYIN, 1 4 )  I NAME  < 1 ) 
forma  t (  S39  ) 

type”The  background  and  image  levels  should  be  an  integer, " 
type"between  0  and  15  .  " 

accepf’ENTER  BACKGROUND  LEVEL  . >",  IBGL 

accept "ENTER  IMAGE  LEVEL  . >".  IVAL 

type"  " 

type" - program  running  - " 

CREATE  BUFFER 

call  pfmt< IBUF,  256,  256,  4,  1  ) 
call  makb(IBUF) 

FILL  BUFFER  WITH  BACKGROUND 
do  20  i=  1,256 

I ARRAY  < a  IDOL  , 
continue 
do  30  i =  1 , 256 

call  prom (I BUF ,  i ,  1,256,  I  ARRAY ) 
continue 


30 


PRIMITIVES  TO  CREATE  THE  IMAGE 

call  C  I RCLE  <  IB'JE,  IVAL-  0,  0,  10) 

call  CIRCLE  < I BIT,  IVAL, O-  0, 20) 

call  CIRCLE < IBUF, IVAL, 0, 0, 25) 
call  Cl RCLE < IDLE. IVAL, 0,  0,  26) 
call  CIRCLE< IBUF, IVAL. 0. O, 27) 

call  LINE  <  IBUF  ,  IVAL,  -1  10,  110,  110,  110) 
call  LINE  (  IBUF ,  IVAL,  -105,  105,  105,  105) 
call  L I NE  < IBUF,  IVAL,  -100,  100,  100,  100) 

call  LINE ( IBUF , IVAL, -1 10, -1 10, 1 10, -1 10) 
call  LINE<  IBUF,  IVAL,  -105,-105,  105,-105) 
call  L 2 NE <  IBUF,  IVAL,  -100,  -100,  100,  -100) 

call  L I  NE  (  IBUF,  IVAL,  1  10,  -1 10,  110,  110) 
call  LINE  (  IBUF,  IVAL,  105,  -105,  105,  105) 
call  L I NE ( IBUF, IVAL, 100, -100, 100, 100) 

call  LINEC IBUF,  IVAL,  -1  10,  -1 10,  -1 10,  110) 
call  LINE  <  IBUF,  IVAL,  -105,  -105,  -105,  105) 
call  LI NEC IBUF, IVAL, -100, -100, -100, 100) 

call  LINEC IBUF, IVAL, -89, 86, 89, 86) 
call  LINEC IBUF,  IVAL,  -86,  83,  86,  83) 
call  LINEC IBUF, IVAL, -83, 80, 83, 80) 

call  LINEC IBUF, IVAL, -89, -86, 89, -86) 
cal  1  LINEC IBUF,  IVAL,  -86, -83, 86,  -83) 
call  LINEC IBUF,  IVAL,  -83,  -80, 83, -BO) 


call  LINEC IBUF,  IVAL,  84.  -71, 84.  71 ) 
call  LINEC  IBUF,  IVAL,  87,  -74,  87,  74) 
call  LINEC IBUF, IVAL, 90, -77, 90, 77) 

call  LINEC IBUF, IVAL, -84, -71, -84, 71 ) 
call  LINEC IBUF, IVAL, -87, -74, -87, 74 ) 
call  LINEC IBUF, IVAL, -90, -77, -90, 77) 

call  LINEC IBUF, IVAL, -50, 50, 50, 50) 
call  LINEC IBUF,  IVAL,  50,  -50,  50,  50) 
call  LINEC IBUF, IVAL, -50, -50, 50, -50) 
call  LINEC IBUF, IVAL, -50, -50, -50, 50) 
call  LINEC  IBUF,  IVAL,  -50,  50,  50,  -50) 
call  LINE  (IBUF,  IVAL,  -50,  -50,  50,  50) 

call  p icout ( IBUF, 1NAME, 1 ) 
call  relb < IBUF) 
s  t  op"  C7X7X7D- image.  fr" 
end 


This  subroutine  is  used  to  draw  circles  in 
spatial  domain  i  Ref.  Chapter  4  ). 

subroutine  circle<IBUF, I VAL, XI , Y1 , IR ) 
integer  IBUF<300>.  IVAL,  XI,  Yl,  IR 
integer  X,  Y,  D.  IROW,  ICOL 

X=  O 
Y=  IR 

D=  3  -  2*IR 

ICOL=  <X+X1)  +  128 
IROW= ( -1 ) * ( Y+Yl )+128 

call  ppnt<IBUF, IROW, ICOL, IVAL) 

ICOL=<Y+Xl >4128 

IR0W=<-1 )*(X4Y1 >4128 

call  ppnt < IBUF, IROW, ICOL, IVAL) 

ICOL- (Y4Y1 >4128 

IROW=  <  — 1 )*< -X4Y1 >4128 

call  ppnt( IBUF, IROW, ICOL, IVAL) 

ICOL=<  X4X1 >4128 

IROW- <  —  1 >*(-Y4Yl >4128 

call  ppnt< IBUF. IROW, ICOL, IVAL) 

ICOL=(-X4Xl >4128 

I ROW= ( — 1 )  *  (  — Y  +Y 1 >4128 

call  ppnt( IBUF,  IROW,  ICOL,  IVAL) 

ICOL=(-Y4Xl >4128 

IR0W=<-1 ) *< -X4Y1 >4l2B 

call  ppnt< IBUF, IROW, ICOL, IVAL) 

ICOL=(-Y4Xl >4128 

IROW=(-l )*(X4Y1 >4128 

call  ppnt< IBUF, IROW, ICOL, IVAL) 

ICOL=  <  — X  +  X 1 >4128 

IROW=(-l >*<Y4Y1 >  4128 

call  ppnt< IBUF, IROW, ICOL, IVAL) 

i  f  (  X.  eq.  Y )  call  p  pnt  (  I BUF,  IROW,  ICOL,  IVAL  > 

if  <  X.  g  e.  Y )  go  to  lb 

if  (D.  ge  0)  go  to  12 
D=  D  +  4*X  +  6 
X=  XM 
go  to  10 

D=  D4  ( X-Y  >  4  10 
Y=Y— 1 

X  =  X4l 

go  to  10 

return 

end 


T 


This  Subroutine  is  used  to  dram  line  in  the 
spatial  domain  (ref  Chapter  4) 


subroutine  LI  NE  (  lbuf-  ivsL  xl<  y  1 .  x  2 .  y2) 

integer  i  b  u  f .  ival,  x  1  >  y  1 .  x  2.  y2 

integer  lx.  ly.  inrrl.  incr2.  sx.  sy.  Id.  xend 

integer  intx.inty.dm 

real  dx.dy.x.y.m 

if(xl.  eq.  x2)  go  to  300 


d  y  =  y 2-y 1 
d  x  =  x2-x 1 
m=d  y /d  x 

i  f  <  m.  1 1.  0 )  go 
1 x  =  x2-x 1 
ly=y2-y 1 
if  <  1  x.  It.  O)  1  x  = 
i  f  <  ly.  It  O)  1  y  = 
1 d=2*l y-1 x 
incr 1=2*1 y 
i ncr2=2* ( 1 y -1 x ) 


to  150 


1 x= ( -1 ) *1 x 
ly= (-1 )*ly 


if  (xl.  It  x2)  go  to  25 
s  x -  x2 
sy  =  y2 
x  end  =  x 1 
go  to  30 
sx  =  x  1 
sy  =  y  1 
x  end  =  x2 

ir  ou»=  (  -1  )  *s  y +  1 28 
icol=sx+128 

call  ppnt(ibuf.  irou.  icol.  ival) 

if  (sx  ge  xend)  go  to  50 
s  x  =  s  x  + 1 

if (Id  ge. 0)  go  to  40 

1  d  =  l d  +  incr 1 

go  to  45 

sy=sy +1 

ld=ld+incr2 

ir  ou»=  (  -1  )  *sy  +  128 

icol=sx+128 

call  ppnt(ibuf.  irou.  icol.  ival) 
go  to  35 


V 


if(dx  lt.dy)  go  to  500 

b  =  q  1  -m<  x  1  I*, 

do  160  i  - x  1  >  x2  ^ 

X  =  1  *1 
y  =m*-  x  +b 

i nty =y 

int  x  =  x 

ir  oui=  < -1  )  *inty  -*-128 
icol=  intx+128 

call  ppnt<ibuf>  irow-  icol>  ival) 


continue 
go  to  50 
x  =  x  1 

d  o  200  i  =  y  1 ,  y 2 
int x  =  x 

irou)=<  -1  )  *i  +  128 
icol=intx+  128 

call  ppnt<ibuf» irow<  icol» ival) 
x  =  x  -M 1/m) 
continue 
go  to  50 

if  (  yl.  It.  y2)  go  to  350 
dm=y  1 
y  l=y2 
y  2=dm 

do  250  i=y 1 ,  y 2 

irow= ( -1 )*i  +  129 
i c  o 1 =  xl  +  128 

call  ppnt(ibuf,  irow.  icol>  ival) 
continue 
return 


This  subroutine  creates  a  subimage  from  an  image 
INAME  :  input  file  name 
I OUTNAME  :  output  file 
ISIZE  :  output  file  size 

C,  R  starting  column  and  row  in  the  input  image 

SUBROUTINE  SUBIMG< INAME, IOUTNM, ISIZE, C, R) 

PARAMETER  NBUFSZ=300 

PARAMETER  NC0LSMAX=256 

INTEGER  I BUF2 ( NBUFSZ  > ,  IBUF(NBUFSZ) 

INTEGER  I ARRAY  <  NCOLSMAX ) 

INTEGER  C,  R 

PICBUF  FOR  INPUT  FILE 

CALL  PICFMT( IBUF, INAME, IHDR ) 

CALL  GFMT < IBUF, NROWS,  NCOLS, NB ITS,  IMODE) 

CALL  MAKB(IBUF) 

CALL  PICIN( IBUF, INAME, IHDR ) 


DELETE  OUTPUT  FILE  IF  IT  EXISTS 
CALL  DFILW( IOUTNM, IER) 

CALL  CHECK (IER) 

I F  (  IER.  NE.  1  )  GO  TO  1  5 
CONTINUE 

I FACTOR- 256/ ISIZE 
PICBUF  FOR  OUTPUT  FILE 

CALL  GFMT  < IBUF,  NROWS,  NCOLS,  NB ITS,  IMODE) 
NROWS2=NROWS/ IFACTOR 
NC0LS2=NC0LS/ IFACTOR 

CALL  PFMT ( IBUF2,  NR0WS2,  NC0LS2.  NB ITS,  IMODE  ) 
CALL  MAKB ( IDUF2 ) 

ISTART=R 
I STOP=R+ I S I ZE - 1 
ICOL=C 
IR0W2=1 
I COL 2= 1 

DO  10  I ROW= I  ST  AR  T ,  ISTOP 

CALL  GROW < I BUF ,  IROW,  ICOL.  ISIZE,  1  ARRAY ) 

CALL  PROW  < I BUF2 ,  IR0W2,  IC0L2,  ISIZE,  I  ARRAY ) 

IR0W2=IR0W2M 

CONTINUE 


CALL  PICOUT< IBUF2,  IOUTNM,  IHDR ) 
CALL  RELB ( I BUF2 ) 

CALL  RELB (IBUF) 

RETURN 


no  ooo  o  o  o  no 


This  subroutine  converts  an  integer  fil 
to  a  complex  form 

ISIZE  :  size  of  the  files 
INNAME,  integer  file  name 
OUTNAME  :  complex  (output)  file  name 


SUBROUTINE  SVTC(ISIZE.  INNAME, OUTNAME ) 
INTEGER  TMP(256>, IN( 1024) 

INTEGER  I NNAME <  7 )  > OUTNAME  <  7 ) 

COMPLEX  OUT (1024) 

CALL  0PEN(S, INNAME, 1, IER > 

CALL  CHECK ( IER ) 

CALL  DFILW ( OUTNAME, IER) 

CALL  CFILW (OUTNAME, 2, KER ) 

CALL  CHECK (KER) 

CALL  0PEN(9, OUTNAME, 2, IER ) 

CALL  CHECK (IER) 

OPEN  9,  OUTNAME, ATT= "OR", ERR= 100 
CONTINUE 

IFACTOR=256/ISIZE 
ITIMES=64/IFACT0R**2~1 

DO  30  1=0, I TIMES 
CALL  RDBLK(8,  I,  TMP,  1,  IER) 

CALL  CHECK (IER) 

CALL  UNP (256, TMP, IN) 

DO  40  J=l, 1024 

OUT ( J ) =CMPLX ( FLOAT ( IN( J) ), O  0) 

40  CONTINUE 

CALL  WRBLK  (  9,  16*1,  OUT,  16.  IER) 

CALL  CHECK < IER) 

30  CONTINUE 

CALL  RESET 
RETURN 
END 
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non 


«r. 


c 

c 

c 

c 

c 

c 

c 

c 

c 


30 


99 


C 

C 

C 

C 


(Ref.  KING,  D.  A.  ,  INHOUSE  ARCHIVE  ) 

(AFIT-  Digital  Signal  Processing  Lab.  ) 

OPTIONS  : 

1  Forward  transform  (Origin  at  CENTER) 

2  Inverse  transform  (Origin  at  CENTER) 

3  :  Forward  transform  (Origin  at  CORNER) 

4  Inverse  transform  (Origin  at  CORNER) 

5  Transform  2D  complex  file 


SUBROUTINE  SFFT2 < I NAME,  OPT  I ON ) 

COMPLEX  X ( 512 ) ,  Y<512> 

INTEGER  OPTION 

I CHAN= 1 

CALL  OPEN< ICHAN, INAME, 2, IER ) 

IF  (  IER.  EQ.  1  )  GO  TO  30 

TYPE "in  SFFT2  err or— op en i n i n  file", IER 
CALL  CLOSE ( ICHAN, IER ) 

GO  TO  99 

COMPUTE  2D-FFT 

I OPT=OPT I ON 

CALL  FFT2UCHAN.  X,  Y,  IOPT,  IER) 

IF  (  IER.  EQ.  1*)  TYPE  "NORMAL  COMPLETION" 
IFIIER  EQ  4)  TYPE " I NVAL I D  FILE  SIZE" 
RETURN 
END 


This  subroutine  will  unpack  four 
from  a  16- bit  integer  word  The 
have  to  be  unpacked  if  each  pi«e 
separately  <  ref  ) 

SUBROUTINE  UNP  < N,  P I XWORD.  PIXELS) 
INTEGER  PIXWORD(N),  PIXELS<4.  N> 

DO  1  I  =  1 ,  N 

DO  1  J= 1 . 4 

P I XELS ( ( 5- J ) , I > - l 5  AND  PIXWORD(I) 
P I XWORD (  I  )  =  I  SHF T ( P I XWORD (  I  )  .  -4  >  , 

I F ( MOD ( N, 12)  EQ  0)  TYPE"UNP  DONE 
RETURN 
END 


4-bit  integers 
lxels  in  a  video  f  i  1  < 
is  to  be  operated  on 


four  pixels  per  wot d 
N  '  a  1  1  ows  highei  order 
arrays  to  be  passed 

■  pick  off  right  piiel 
shift  word  4  bite  right 

■  N  ,  to  pick  off  fie  x  t  pixel 
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A  line  estimation  method  is  described  to  estimate  the  degrading  OTF  of  the 
linear  camera  motion  in  one  direction.  An  additional  method  called  the 
"logarithmic  estimate"  is  also  developed.  Frequency  and  space 
domain  estimation  methods  are  developed  for  the  noise-to-signal 
ratio  for  the  case  of  additive  gaussian  noise.  Inverse,  Wiener,  and  power 
spectrum  equalization  filters  based  on  these  estimates  are  implemented 
on  a  digital  computer  to  restore  a  variety  of  degraded  images  and  examples 
of  blind  restoration  of  these  degraded  images  are  presented. 
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